In [1]:
# Parameters
quantile = 0.316666666667
SWARM_PATH = "/newhome/kamran.karimi1/Project/git/Swarm/dataset/RidgeCrest/DataS1_noXYZ.txt"

fetch the data: Shelly 2020

In [2]:
#--- extract the catalog (new)

import pandas as pd
#import matplotlib; matplotlib.use('agg')
import sys
from matplotlib import rc
import matplotlib
import matplotlib.pylab as plt
import numpy as np
from matplotlib import patches
import warnings
import matplotlib.ticker
import matplotlib as mpl
#from mpl_toolkits.mplot3d import Axes3D
import datetime
#import path
from math import *
import matplotlib.dates as mdates
from matplotlib import dates
import matplotlib.dates as md
import matplotlib.cm as cm
from matplotlib.font_manager import FontProperties
import itertools
import traceback
from scipy.ndimage import gaussian_filter
import math

def DrawFrame(ax, (alpha_xl,alpha_xr),(alpha_yb,alpha_yt),linewidth,LOG_X=None,LOG_Y=None):
    [xlo,xhi,ylo,yhi] = ax.axis()
    if LOG_X:
        [xlo,xhi,junk,junk] = np.log10(ax.axis())
    if LOG_Y:
        [junk,junk,ylo,yhi] = np.log10(ax.axis())
    lx = xhi - xlo
    ly = yhi - ylo
    xy = [xlo - alpha_xl * lx, ylo - alpha_yb * ly]
    height = ly*(1+alpha_yb+alpha_yt)
    width = lx*(1+alpha_xl+alpha_xr)
    xy_end=[xy[0]+width,xy[1]+height]
    if LOG_X:
        xy[0] = 10 ** xy[0]
        xy_end[0] = 10 ** xy_end[0]
    if LOG_Y:
        xy[1] = 10 ** xy[1]
        xy_end[1] = 10 ** xy_end[1]
    ax.add_patch( patches.Rectangle(xy=xy, width=xy_end[0]-xy[0], 
                                    height=xy_end[1]-xy[1], linewidth=linewidth,
                                    clip_on=False,facecolor=None,edgecolor='black',fill=None) ) 
    
#--- add a new time attribute
def ConvertTime( df_in ):
    df=df_in.copy()
    df.insert(0,'date',pd.to_datetime(swarm[['year', 'month', 'day', 'hour', 'minute', 'second']]))
    df.drop(['year', 'month', 'day', 'hour', 'minute', 'second'],axis=1,inplace=True)
    return df

#--- set path
#SWARM_PATH = './dataset/RidgeCrest/DataS1_noXYZ.txt' #--- comment if passed by argument
DIR_OUTPT = '.' #'/dataset/RidgeCrest' #sys.argv[2]
DIR_OUTPT_figs = '.' #'/dataset/figs' #'/Users/Home/Desktop/Tmp/txt' #'/Users/Home/Dropbox/Documents/papers/microSeismicityOklahoma/figs'

#--- store
swarm = pd.read_csv( SWARM_PATH, sep = ',' ) #--- parse data

swarm = ConvertTime( swarm ) #--- add new column 'date'

#--- sort based on time
swarm.sort_values(by=['date'],inplace=True)

#--- reindex
swarm.reset_index(inplace=True,drop=True)

#--- rename cols
swarm = swarm.rename(index=str,columns={'lat':'latitude','lon':'longitude','mag':'magnitude'})

swarm.head()
Out[2]:
date latitude longitude depth magnitude ID
0 2019-07-04 15:35:29.400 35.708036 -117.499365 12.527 0.17 70456129
1 2019-07-04 15:42:47.900 35.708040 -117.499382 12.508 0.59 70456568
2 2019-07-04 16:07:19.950 35.708081 -117.499430 12.517 0.50 70458040
3 2019-07-04 16:13:11.390 35.707955 -117.499211 12.506 0.04 70458391
4 2019-07-04 16:13:43.160 35.708044 -117.499382 12.495 1.50 70458423
In [3]:
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings

rc('text', usetex=True)
font = {'size'   : 20}
matplotlib.rc('font', **font)
In [4]:
    
# #--- add a new time attribute
# def ConvertTime( df_in ):
    
#     df=df_in.copy()
#     df['date']=pd.to_datetime(df['Event ID'] +' ' + df['Origin Time'])
#     df['Event ID']=df.index
# #    df.insert(0,'date',pd.to_datetime(swarm[['year', 'month', 'day', 'hour', 'minute', 'second']]))
#     df.drop(['Origin Time'],axis=1,inplace=True)
#     return df

# #--- set path
# SWARM_PATH = './dataset/RidgeCrest/Dataset S2.txt' #--- comment if passed by command line
# DIR_OUTPT = '.' #'./dataset/RidgeCrest' #sys.argv[2]
# DIR_OUTPT_figs = '.' #'./dataset/figs' #'/Users/Home/Desktop/Tmp/txt' #'/Users/Home/Dropbox/Documents/papers/microSeismicityOklahoma/figs'

# #--- store
# swarm = pd.read_csv( SWARM_PATH, sep = ',' ) #--- parse data

# swarm = ConvertTime( swarm ) #--- add new column 'date'

# #--- sort based on time
# swarm.sort_values(by=['date'],inplace=True)

# #--- reindex
# swarm.reset_index(inplace=True,drop=True)

# #--- rename cols
# swarm = swarm.rename(index=str,columns={'Latitude':'latitude','Longitude':'longitude',
#                                         'Magnitude (SCEDC preferred)':'magnitude','Depth (km)':'depth'})

# swarm.head()
In [5]:
#--- remove nan

#swarm = swarm.dropna(subset=['magnitude'])
swarm.info()
<class 'pandas.core.frame.DataFrame'>
Index: 34091 entries, 0 to 34090
Data columns (total 6 columns):
date         34091 non-null datetime64[ns]
latitude     34091 non-null float64
longitude    34091 non-null float64
depth        34091 non-null float64
magnitude    34091 non-null float64
ID           34091 non-null int64
dtypes: datetime64[ns](1), float64(4), int64(1)
memory usage: 1.8+ MB

Moment tensors scedc

In [6]:
    
# #--- add a new time attribute
# def ConvertTime( df_in ):
    
#     df=df_in.copy()    
#     df['date']=pd.to_datetime(df['Date']+' '+df['Time(UTC)'])
# #    df['date']=pd.to_datetime(df['Event ID'] +' ' + df['Origin Time'])
# #    df['Event ID']=df.index
# #    df.insert(0,'date',pd.to_datetime(df_moment[['year', 'month', 'day', 'hour', 'minute', 'second']]))
#     df.drop(['Time(UTC)','(Ml)','Mp','originalDepth','Solutions','VR','Date'],axis=1,inplace=True)
#     return df

# #--- set path
# SWARM_PATH2 = './dataset/RidgeCrest/momentTensor.csv' #--- comment if passed by command line

# #--- store
# df_moment = pd.read_csv( SWARM_PATH2, sep = ',' ) #--- parse data

# df_moment = ConvertTime( df_moment ) #--- add new column 'date'

# #--- sort based on time
# df_moment.sort_values(by=['date'],inplace=True)

# #--- reindex
# df_moment.reset_index(inplace=True,drop=True)

# #--- rename cols
# df_moment = df_moment.rename(index=str,columns={'EventID':'Event ID','Lat.':'latitude',
#                                                 'Long.':'longitude','Mw':'magnitude'})

# df_moment.head()

fetch tensors from the html files

  • use event id's from the above catalog
In [7]:
# azimuth_t=[]
# azimuth_n=[]
# azimuth_p=[]
# plunge_t = []
# plunge_n = []
# plunge_p = []
# mt = []
# mn = []
# mp = []
# for items in df_moment.itertuples():
#     if 1: #try:
#         eventID = items[1]
#         html_url='https://service.scedc.caltech.edu/MomentTensor/solutions/web_%s/ci%s_MT.html'%(eventID,eventID)
# #        print html_url
#         dff=pd.read_html(html_url)
#         dff=dff[2] #--- eigenvalues and eigenmodes
        
# #        if eventID not in list(swarm['Event ID']):
# #            continue
# #         #--- stress catalog
# # #        if eventID == 38636511:
# #         print eventID, df_st[df_st['Event ID']==eventID]['Origin Time'].iloc[0]-items[4]
# #            print df_st[df_st['Event ID']==eventID]['Origin Time'].iloc[0],items[4]

# #         #--- parse moment tensor
#         azimuth_t.append(float(dff.iloc[2][3]))
#         azimuth_n.append(float(dff.iloc[3][3]))
#         azimuth_p.append(float(dff.iloc[4][3]))
#         plunge_t.append(float(dff.iloc[2][2]))
#         plunge_n.append(float(dff.iloc[3][2]))
#         plunge_p.append(float(dff.iloc[4][2]))
    
    
#         mt.append(float(dff.iloc[2][1]))
#         mn.append(float(dff.iloc[3][1]))
#         mp.append(float(dff.iloc[4][1]))
#         assert float(dff.iloc[2][1]) > 0.0
#         assert float(dff.iloc[4][1]) < 0.0

insert as columns

In [8]:
# df_moment.insert(len(df_moment.keys()),'mt',mt)
# df_moment.insert(len(df_moment.keys()),'Azimuth(T)',azimuth_t)
# df_moment.insert(len(df_moment.keys()),'Plunge(T)',plunge_t)

# df_moment.insert(len(df_moment.keys()),'mn',mn)
# df_moment.insert(len(df_moment.keys()),'Azimuth(N)',azimuth_n)
# df_moment.insert(len(df_moment.keys()),'Plunge(N)',plunge_n)

# df_moment.insert(len(df_moment.keys()),'mp',mp)
# df_moment.insert(len(df_moment.keys()),'Azimuth(P)',azimuth_p)
# df_moment.insert(len(df_moment.keys()),'Plunge(P)',plunge_p)

common events between the two catalogs

In [9]:
# df_moment.index=df_moment['Event ID']
# common_elm = list(set(df_moment['Event ID']) & set(swarm['Event ID']))
# df_moment = df_moment.loc[common_elm]
# df_moment.reset_index(drop=True, inplace=True)
# df_moment.head()

moment tensors scedc

Column legend:

EVID : Event Id
ET : Earthquake Type
MAG : Event Magnitude
M : Magnitude Type
LAT : Latitude
LON : Longitude
DEPTH : Earthquake Depth
Q : Quality
NPH : Number of phases or arrivals
WRMS : Weighted RMS
ERHOR : Horizontal Error (location uncertainty that accompanies the location)
ERDEP : Depth Error (standard deviation of depth estimate)
ERTIME : Origin Time Error (standard deviation of origin time)
STRIKE : Strike of fault plane
DIP : Dip of fault plane
RAKE : Rake of fault plane
FPUC : Fault Plane Uncertainty in degrees
APUC : Auxialiary Plane Uncertainty in degrees
NPPL : Number of P Polarities
MFRAC : Fraction Misfit Polarities
FMQ : Focal Mechanism Quality
PROB : Probability true mechanism is "close" to preferred solution
STDR : Station Distribution
NSPR : Number of S/P Ratios
MAVG : Average log10(S/P) Misfit

In [10]:
    
# #--- add a new time attribute
# def ConvertTime( df_in ):
    
#     df=df_in.copy()
# #    display(df.head())
#     df.insert(loc=2, column='date', value=pd.to_datetime(df['#YYY/MM/DD']+' '+df['HH:mm:SS.ss']))
# #    df['date']=pd.to_datetime(df['Event ID'] +' ' + df['Origin Time'])
# #    df['Event ID']=df.index
# #    df.insert(0,'date',pd.to_datetime(df_moment[['year', 'month', 'day', 'hour', 'minute', 'second']]))
#     df.drop(['#YYY/MM/DD','HH:mm:SS.ss', u'ET', u'GT',u'M', u'Q', u'NPH', u'WRMS', u'ERHOR', u'ERDEP',
#        u'ERTIME', u'STRIKE', u'DIP', u'RAKE', u'FPUC', u'APUC', u'NPPL',
#        u'MFRAC', u'FMQ', u'PROB', u'STDR', u'NSPR', u'MAVG'],axis=1,inplace=True)
#     return df

# #--- set path
# SWARM_PATH3 = './dataset/RidgeCrest/catalog_fmsearch.csv' #--- comment if passed by command line

# #--- store
# df_moment = pd.read_csv( SWARM_PATH3, sep = ',' ) #--- parse data

# df_moment = ConvertTime( df_moment ) #--- add new column 'date'

# #--- sort based on time
# df_moment.sort_values(by=['date'],inplace=True)

# #--- reindex
# df_moment.reset_index(inplace=True,drop=True)

# #--- rename cols
# df_moment = df_moment.rename(index=str,columns={'EVID':'Event ID','LAT':'latitude',
#                                                 'LON':'longitude','MAG':'magnitude', 'DEPTH':'depth'})

# #df_moment = df_moment.iloc[:32]
# df_moment.head()
# #df_moment.keys()
In [11]:
#df_moment.info()

fetch angles

In [12]:
# def func(x):
#     try:
#         html_url = x
#         dff=pd.read_html(html_url)
#         dff=dff[3] #--- eigenvalues and eigenmodes
#         #--- parse moment tensor
#         assert dff.iloc[0]['Axis'] == 'P'
#         assert dff.iloc[1]['Axis'] == 'T'
#         assert dff.iloc[1]['Azimuth'][-4:-1] == '&de'
#         assert dff.iloc[0]['Azimuth'][-4:-1] == '&de'
#         assert dff.iloc[1]['Plunge'][-4:-1] == '&de'
#         assert dff.iloc[0]['Plunge'][-4:-1] == '&de'
#         azimuth_t = float(dff.iloc[1]['Azimuth'][:-4])
#         azimuth_p = float(dff.iloc[0]['Azimuth'][:-4])
#         plunge_t  = float(dff.iloc[1]['Plunge'][:-4])
#         plunge_p  = float(dff.iloc[0]['Plunge'][:-4])
#         return (azimuth_t,plunge_t,azimuth_p,plunge_p)        
#     except:
#         return (None, None, None, None )
# df_moment['html']=df_moment['Event ID'].apply(lambda x: 'https://service.scedc.caltech.edu/FocMech/ci%s.cifm1.html'%(x))
# s_series = df_moment['html'].apply(func)
# df_moment['Azimuth(T)'] = s_series.apply(lambda x: x[0])
# df_moment['Plunge(T)'] = s_series.apply(lambda x: x[1])
# df_moment['Azimuth(P)'] = s_series.apply(lambda x: x[2])
# df_moment['Plunge(P)'] = s_series.apply(lambda x: x[3])

# df_moment.drop(['html'],axis=1,inplace=True)
In [13]:
#df_moment.head() 

save data frame

In [14]:
#df_moment.to_pickle('./dataset/RidgeCrest/momentTensors.pkl')
#swarm=pd.read_pickle('./dataset/RidgeCrest/momentTensors.pkl')
In [15]:
#--- plot spatial map 
#--- fault map california: temporal evolution of events


def DataFrameSubSet( df, column, (xmin,xmax) ):
    return df[ ( df[column] >= xmin ) & 
               ( df[column] < xmax ) ]

#--- subset of data
swarm_lohi = swarm.copy()#DataFrameSubSet( swarm, 
                          #   'date', 
                           #  ( pd.to_datetime('2010-04-04'), pd.to_datetime('2010-06-26') ) )

swarm_lohi.plot.scatter('longitude','latitude',
                        s=3**(swarm_lohi['magnitude']),
                        c='date',cmap='jet',
                        alpha=0.1,figsize=(5,4)) #--- plot

# swarm_lohi.plot.scatter('longitude','depth',
#                         s=3**(swarm_lohi['magnitude']),
#                         c='date',cmap='jet',
#                         alpha=0.1,figsize=(5,4)) #--- plot

# swarm_lohi.plot.scatter('latitude','depth',
#                         s=3**(swarm_lohi['magnitude']),
#                         c='date',cmap='jet',
#                         alpha=0.1,figsize=(5,4)) #--- plot

plt.show()
#plt.savefig('%s/map.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')

#plt.figure(figsize=(6,6)).gca(projection='3d')
#plt.xlabel('Long')
#plt.ylabel('Lat')
#plt.scatter(swarm_lohi['longitude'],
 #           swarm_lohi['latitude'],
  #          swarm_lohi['depth']) #s=10*swarm['magnitude'],cmap='jet') #--- plot

make movies

In [16]:
# import os

# class makeDATA:
#     def __init__(self, write, read = None ):
#         self.wirite = write
#         self.read = read   # embedded object: attrs (coord, lohi, type, vel)

#     def __del__( self ):
#         self.wirite.close()

#     def Map(self, coord, comnt, *args ):
#          nAtom = len(coord)
#          self.wirite.write('%s\n%s\n'%(nAtom,comnt))
#          for i in coord:
#              ndime = len( coord[i] )
#              self.wirite.write('%s\t'% i)
#  #           self.wirite.write('%d\t'% type[i])
#              for idime in xrange( ndime ):
#                  self.wirite.write('%e\t'%coord[i][idime] )
#              if args:
#                  for dics in args:
#  #                   if type( dics[ i ] == type( 1 ) ):
#  #                       self.wirite.write( '%d\t' % dics[ i ] )
#  #                   else:
#                      VAL = 0.0 if abs(dics[ i ])<1.0e-15 else dics[ i ]
#                      self.wirite.write( '%s\t' % VAL )
#                  self.wirite.write( '\n' )
#              else:
#                  self.wirite.write( '\n' )
#          self.wirite.flush()
#          self.wirite.close()

# swarm_lohi.sort_values(by='date',inplace=True)            
# coord_length=swarm_lohi['longitude'].apply(lambda x:[x])+swarm_lohi['latitude'].apply(lambda x:[x])
# Magnitude = np.exp(swarm_lohi['magnitude'])/20000

# makeDATA( open( '%s/map.xyz'%DIR_OUTPT, 'w' ), None ).Map( coord_length.to_dict(), 'ITIME=%s'%0, Magnitude.to_dict()) #--- instantanuous
# EVENT_ID = Magnitude.keys()
# EVENT_ID = map(int,EVENT_ID)
# EVENT_ID.sort()
# cords={}
# mag={}
# t0 = None
# os.system('rm %s/temporalMap.xyz'%DIR_OUTPT)
# for j_event,indexx in zip(EVENT_ID,xrange(10000)):
# #    print j_event
#     t = swarm_lohi['date'].iloc[ j_event ]
#     if not t0:
#         t0 = t
#     assert t >= t0, 't0=%s,t=%s'%( t0, t )
#     #--- temporal map 
#     cords[ j_event ] = coord_length[ j_event ]
#     mag[ j_event ] = Magnitude[ j_event ]
#     itime=swarm_lohi[ 'date' ].iloc[ j_event ]
#     tstr='%s-%s'%(itime.strftime("%x"), itime.strftime("%X"))
#     makeDATA( open( '%s/temporalMap.xyz'%DIR_OUTPT, 'a' ), None ).Map( cords, 'ITIME=%s'%itime, mag ) #--- instantanuous
#     t0 = t
    
    
    
In [17]:
# #--- plot timeseries


class Tlohi:
     def __init__(self,lo,hi):
        self.lo = lo
        self.hi = hi

# fig=plt.figure(figsize=(8,4))
# ax=fig.add_subplot(111)
# plt.xlabel('Time')
# plt.ylabel('M')

# plt.ylim(swarm['magnitude'].min(),swarm['magnitude'].max())

# tt=swarm[swarm['magnitude']==swarm['magnitude'].sort_values(ascending=False).iloc[3]]['date']
# plt.xlim(tt-datetime.timedelta(days=30),tt+datetime.timedelta(days=30))# pd.to_datetime('2014-12-04'),pd.to_datetime('2015-06-1'))
# #plt.xlim(swarm['date'].min(),pd.to_datetime('2010-07-20'))
# plt.scatter(swarm['date'],swarm['magnitude'],
#             s=2*np.exp(swarm['magnitude']),
#             alpha=0.2)

# #myTlohi = Tlohi( pd.to_datetime('2010-04-08'), 
# #                 pd.to_datetime('2010-07-01'))

# #plt.plot([myTlohi.lo,myTlohi.lo],
# #        [-2,7],'r-')

# #plt.plot([myTlohi.hi,myTlohi.hi],
# #        [-2,7],'r-')
# fig.autofmt_xdate()
# plt.savefig('%s/timeSeries.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
# plt.show()
In [18]:
#--- freq-magnitude relation

import sys 

def histogramACCUMLTD( slist ):
    assert type( slist ) == type( [] ), 'arg must be a list. a %s is given!' %( type( slist ) )
    d = {}
    for item in slist:
        try:
            d[ item ] += 1
        except:
            d[ item ] = 1
    keys = d.keys()
    keys.sort()

    cdf = 0.0
    xi = min( slist ) - 1.0e-6
    xf = max( slist ) + 1.0e-6
    npoin = len( slist )
    adict = {}
    for ikey, index in zip( keys, xrange( sys.maxint ) ):
        adict[ index ] = [ xi, ikey, cdf ]
        cdf += 1.0 * d[ ikey ] # / npoin
        xi = ikey
    adict[ index + 1 ] = [ xi, xf, cdf ]
    return adict





#--- set min/max time to avoid temporal incompletenesss issue
swarm_copy = swarm.copy()#DataFrameSubSet( swarm, 
                          #   'date', 
                           #  ( myTlohi.lo, 
                            #   myTlohi.hi ) )
#--- accumulated histogram
N = len(swarm_copy['magnitude'])
slist=np.array(swarm_copy['magnitude'])
slist.sort()
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()

#--- plot
fig= plt.figure(figsize=(3,3))#,dpi=150)
ax = fig.add_subplot(111)

ax.axis([-1,6,.9,1e5])
ax.set_yscale('log')

#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(ax.axis()[0],ax.axis()[1],int(nbin)+1))

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))

#--- put minor bins
#locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
#ax.xaxis.set_minor_locator(locmin)
#ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')


xmin=swarm_copy['magnitude'].min()
xmax=swarm_copy['magnitude'].max()
dx=0.1
junk = ax.hist( swarm_copy['magnitude'],
                bins=int((xmax-xmin)/dx),
               label='histogram',color='silver') #--- histogram
xx=[];yy=[]
for ikey in keys:
    xx.append(d[ikey][0])
    yy.append(N-d[ikey][2])
ax.plot(xx,yy,
            linestyle='-', drawstyle='steps-post',color='black',
             linewidth=1.0) #--- accumulated
    
b=0.8
ax.plot([2.5,6],[1e4, 1e4* 10**(-b*(6-2.5))],'r-.',linewidth=1)

DrawFrame(ax, (0.2,0.06),(0.32,0.1),0.01,LOG_Y=True) 

plt.savefig('%s/gr.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')#,transparent=True)
plt.show()
In [19]:
#--- estimate b by varying the mc value

from math import *
import scipy.stats

#--- set min/max time to avoid temporal incompletenesss issue
swarm_copy = swarm.copy()#DataFrameSubSet( swarm.copy()#, 
                          #   'date', #
                           #  ( myTlohi.#lo, 
                            #   myTlohi.hi# ) )

#--- plot distributions (Nm)
fig = plt.figure(figsize=(12,4))
ax=fig.add_subplot(131)
ax.set_xlabel('$m$')
ax.set_ylabel('$N(mag.>m)$')
xmin = np.floor(swarm_copy['magnitude'].min())
xmax = np.ceil(swarm_copy['magnitude'].max())
ax.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax.set_xlim(xmin,xmax)
#plt.ylim(1,1e5)
ax.set_yscale('log')
for ikey in keys:
    ax.plot([d[ikey][0],d[ikey][1]],
             [N-d[ikey][2],N-d[ikey][2]],
             '-',color='black') #--- accumulated
junk = ax.hist( swarm_copy['magnitude'],
                bins=64,
               label='histogram') #--- histogram

#--- list of magnitudes
mm_list = swarm_copy['magnitude'].to_list()
mm_list = np.sort(mm_list)

#--- b vs mc
ax2=fig.add_subplot(132)
ax2.set_xticks(np.linspace(-1.0,3.0,5))
#ax2.set_ylim(0,1.8)
#ax4 = ax2.twinx()
#ax4.set_yscale('log')
#ax4.set_ylim(1e-16,)
#plt.xlim(0.0,1.5)
ax2.set_xlabel('$m_c$')
ax2.set_ylabel('$b$')
#ax4.set_ylabel(r'$p_{val}$',labelpad=-70)

#--- bvalue bhatacharya eq. 5
DM = 0.1
n = 30
#--- 
MC = []
B = []
Err_B = []
pval = []
for mc_new in [0.05+i*0.1 for i in xrange(n)]: #--- m values are quantized (+- 0.1)!
    mlist = swarm_copy[ swarm_copy['magnitude']-mc_new>=0.0]['magnitude'] 
    m_mean=mlist.mean()
    assert m_mean-mc_new > 0.0, 'm_mean=%s,mc=%s'%( m_mean, mc )
    k_mean = (m_mean-mc_new)/DM
    b_value = log(np.exp(1),10)*log(1+DM/(m_mean-mc_new))/DM
    nsamp = len(mlist)
#    err_b = 1.96*(b_value/nsamp/log(10,exp(1.0)))**0.5
    err_b = 2*(k_mean*(k_mean+1)/nsamp)**0.5/log(10,exp(1.0)) #--- Guttorp and Hopkins (1986) eq.(4)

    
    
#     #--- bvalue bhatacharya eq. 6
#     #--- accumulated histogram
#     slist = mm_list[ mm_list>=mc_new ]
#     slist.sort()
#     NN = len(slist)
#     cdf = histogramACCUMLTD( slist.tolist() )
#     #--- bvalue bhatacharya eq. 6
#     n_cdf=len(cdf)-1
#     bi = []
#     for i in xrange(n_cdf):
#         mi = cdf[ i ][ 0 ]
#         Ni = NN - cdf[ i ][ 2 ]
#         bij = []
#         for j in xrange(n_cdf):
#             mj = cdf[ j ][ 0 ]
#             Nj = NN - cdf[ j ][ 2 ]
#             if i == j: continue
#             slope = ( log( Nj, 10 ) - log( Ni, 10 ) ) / ( mj - mi )
#             assert not np.isnan(slope)
#             bij.append( slope )
#         bi.append( np.median( bij ) )
#     bval_2nd = -np.median( bi )
#     err_b2nd = np.std(bi)/(len(bi))**0.5
    
    #--- plot b vs mc
    ax2.errorbar([mc_new],[b_value],yerr=err_b,marker='o',color='black')
#    ax2.errorbar([mc_new],[bval_2nd],yerr=err_b2nd, marker='o',color='red')
    
    #---
    MC.append(mc_new)
    B.append(b_value)
    Err_B.append(err_b)

    #--- p-value (by hand)
#    (dmax,pval2nd) = scipy.stats.kstest( mlist, lambda x: 1.0-10**(-b_value*(x-mc_new)))
#    ax4.plot([mc_new],[pval2nd],marker='o',color='red')
#    pval.append(pval2nd)

#--- choose the range of mc
mc_range=[1.5,1.5]
index_tf = [i and j for i, j in zip(mc_range[0]<=np.array(MC),np.array(MC)<mc_range[1])]
bvall =  np.array(B)[index_tf]
Err_B =  np.array(Err_B)[index_tf]

#--- b vs mc: vertical lines
#ax2.plot([mc_range[0],mc_range[0]],[np.min(B),np.max(B)],'--r')    
#ax2.plot([mc_range[1],mc_range[1]],[np.min(B),np.max(B)],'--r')    

#--- N vs m: vertical line
#ax.plot([mc_range[0]+DM/2,mc_range[0]+DM/2],
 #       [1,len(swarm_copy['magnitude'])],'--r')    

#--- n*10*bm (rescaled)
bval=0.8 #np.mean(bvall)
#bval2=0.8
ax3=fig.add_subplot(133)
ax3.set_xlabel('$m$')
#ax3.set_ylabel(r'$N\times10^{%s\times m}$'%bval)
ax3.set_ylabel(r'$N\times10^{%s\times m}$'%bval)
ax3.yaxis.set_label_position("right")
ax3.yaxis.tick_right()
ax3.set_yscale('log')
ax3.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax3.set_xlim(xmin,xmax)
#ax3.set_ylim(1e3,1e6)
#ax3.set_xticks(np.linspace(ax3.axis()[0],ax3.axis()[1],4))
for ikey in keys: #--- rescaled
    c = 10**(bval*d[ikey][0])
    ax3.plot([d[ikey][0],d[ikey][1]],
             [c*(N-d[ikey][2]),c*(N-d[ikey][2])],'-',color='black',
            ) #--- accumulated
#     c = 10**(bval2*d[ikey][0])
#     ax3.plot([d[ikey][0],d[ikey][1]],
#              [c*(N-d[ikey][2]),c*(N-d[ikey][2])],'-',color='red',
#              ) #--- accumulated
ax3.plot([mc_range[0]+DM/2,mc_range[0]+DM/2], #--- vertical line
        [len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])/10,
         len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])],'--r')    
#plt.plot(MC, pval,'-o',color='black')

fig.tight_layout(w_pad=-0.5)
#fig.savefig('%s/b.png'%DIR_OUTPT_figs,dpi=75)
plt.show()

AttributeErrorTraceback (most recent call last)
<ipython-input-19-0764295145f9> in <module>()
     30 
     31 #--- list of magnitudes
---> 32 mm_list = swarm_copy['magnitude'].to_list()
     33 mm_list = np.sort(mm_list)
     34 

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/pandas/core/generic.pyc in __getattr__(self, name)
   3612             if name in self._info_axis:
   3613                 return self[name]
-> 3614             return object.__getattribute__(self, name)
   3615 
   3616     def __setattr__(self, name, value):

AttributeError: 'Series' object has no attribute 'to_list'
In [20]:
#--- estimate b by varying the mc value

from math import *
import scipy.stats

#--- set min/max time to avoid temporal incompletenesss issue
swarm_copy = swarm.copy()#DataFrameSubSet( swarm.copy()#, 
                          #   'date', #
                           #  ( myTlohi.#lo, 
                            #   myTlohi.hi# ) )

#--- plot distributions (Nm)
fig = plt.figure(figsize=(12,4))
ax=fig.add_subplot(131)
ax.set_xlabel('$m$')
ax.set_ylabel('$N(mag.>m)$')
xmin = np.floor(swarm_copy['magnitude'].min())
xmax = np.ceil(swarm_copy['magnitude'].max())
ax.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax.set_xlim(xmin,xmax)
#plt.ylim(1,1e5)
ax.set_yscale('log')
for ikey in keys:
    ax.plot([d[ikey][0],d[ikey][1]],
             [N-d[ikey][2],N-d[ikey][2]],
             '-',color='black') #--- accumulated
junk = ax.hist( swarm_copy['magnitude'],
                bins=64,
               label='histogram') #--- histogram

#--- list of magnitudes
mm_list = swarm_copy['magnitude'].to_list()
mm_list = np.sort(mm_list)

#--- b vs mc
ax2=fig.add_subplot(132)
ax2.set_xticks(np.linspace(0,1,5,endpoint=True))
#ax2.set_ylim(0.5,1.0)
ax4 = ax2.twinx()
ax4.set_yscale('log')
ax4.set_ylim(1e-16,)
ax2.set_xlabel('$m_c$')
ax2.set_ylabel('$b$')
ax4.set_ylabel(r'$p_{val}$',labelpad=-10)
ax4.set_xticks(np.linspace(0,1,5,endpoint=True))

#--- bvalue bhatacharya eq. 5
DM = 0.01
n = 200
#--- 
MC = []
B = []
Err_B = []
pval = []
for mc_new in [i*DM for i in xrange(n)]: #--- m values are quantized (+- 0.1)!
    mlist = swarm_copy[ swarm_copy['magnitude']-mc_new>=0.0]['magnitude'] 
    m_mean=mlist.mean()
    assert m_mean-mc_new > 0.0#, 'm_mean=%s,mc=%s'%( m_mean, mc )
    k_mean = (m_mean-mc_new)/DM
    b_value = log(np.exp(1),10)*log(1+DM/(m_mean-mc_new))/DM
    nsamp = len(mlist)
    err_b = 1.96 * b_value / nsamp ** 0.5 #1.96*(b_value/nsamp/log(10,exp(1.0)))**0.5
#    err_b = 2*(k_mean*(k_mean+1)/nsamp)**0.5/log(10,exp(1.0)) #--- Guttorp and Hopkins (1986) eq.(4)

    
    
#     #--- bvalue bhatacharya eq. 6
#     #--- accumulated histogram
#     slist = mm_list[ mm_list>=mc_new ]
#     slist.sort()
#     NN = len(slist)
#     cdf = histogramACCUMLTD( slist.tolist() )
#     #--- bvalue bhatacharya eq. 6
#     n_cdf=len(cdf)-1
#     bi = []
#     for i in xrange(n_cdf):
#         mi = cdf[ i ][ 0 ]
#         Ni = NN - cdf[ i ][ 2 ]
#         bij = []
#         for j in xrange(n_cdf):
#             mj = cdf[ j ][ 0 ]
#             Nj = NN - cdf[ j ][ 2 ]
#             if i == j: continue
#             slope = ( log( Nj, 10 ) - log( Ni, 10 ) ) / ( mj - mi )
#             assert not np.isnan(slope)
#             bij.append( slope )
#         bi.append( np.median( bij ) )
#     bval_2nd = -np.median( bi )
#     err_b2nd = np.std(bi)/(len(bi))**0.5
    
    #--- plot b vs mc
    ax2.errorbar([mc_new],[b_value],yerr=err_b,marker='o',color='black')
#    ax2.errorbar([mc_new],[bval_2nd],yerr=err_b2nd, marker='s',color='green')
    
    #---
    MC.append(mc_new)
    B.append(b_value)
    Err_B.append(err_b)

    #--- p-value (by hand)
    (dmax,pval2nd) = scipy.stats.kstest( mlist, lambda x: 1.0-10**(-b_value*(x-mc_new)))
    ax4.plot([mc_new],[pval2nd],marker='o',color='red')
    pval.append(pval2nd)
#    print b_value, pval2nd 

xmin = min(MC)
xmax = max(MC)
ax2.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/.2))
ax2.set_xlim(xmin,xmax)


#--- choose the range of mc
mc_range=[1.5,1.5]
index_tf = [i and j for i, j in zip(mc_range[0]<=np.array(MC),np.array(MC)<mc_range[1])]
bvall =  np.array(B)[index_tf]
Err_B =  np.array(Err_B)[index_tf]

#--- b vs mc: vertical lines
#ax2.plot([mc_range[0],mc_range[0]],[np.min(B),np.max(B)],'--r')    
#ax2.plot([mc_range[1],mc_range[1]],[np.min(B),np.max(B)],'--r')    

#--- N vs m: vertical line
#ax.plot([mc_range[0]+DM/2,mc_range[0]+DM/2],
 #       [1,len(swarm_copy['magnitude'])],'--r')    

#--- n*10*bm (rescaled)
bval=0.8 #np.mean(bvall)
ax2.plot([xmin,xmax],[bval,bval],'--r')
#bval2=0.8
ax3=fig.add_subplot(133)
ax3.set_xlabel('$m$')
#ax3.set_ylabel(r'$N\times10^{%s\times m}$'%bval)
ax3.set_ylabel(r'$N\times10^{%s\times m}$'%bval)
ax3.yaxis.set_label_position("right")
ax3.yaxis.tick_right()
ax3.set_yscale('log')
xmin = np.floor(swarm_copy['magnitude'].min())
xmax = np.ceil(swarm_copy['magnitude'].max())
ax3.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax3.set_xlim(xmin,xmax)
#ax3.set_ylim(1e3,1e6)
#ax3.set_xticks(np.linspace(ax3.axis()[0],ax3.axis()[1],4))
for ikey in keys: #--- rescaled
    c = 10**(bval*d[ikey][0])
    ax3.plot([d[ikey][0],d[ikey][1]],
             [c*(N-d[ikey][2]),c*(N-d[ikey][2])],'-',color='black',
            ) #--- accumulated
#     c = 10**(bval2*d[ikey][0])
#     ax3.plot([d[ikey][0],d[ikey][1]],
#              [c*(N-d[ikey][2]),c*(N-d[ikey][2])],'-',color='red',
#              ) #--- accumulated
ax3.plot([mc_range[0]+DM/2,mc_range[0]+DM/2], #--- vertical line
        [len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])/10,
         len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])],'--r')    
#plt.plot(MC, pval,'-o',color='black')

# mc_new = 0.5
# b_value = 0.94
# ax2.scatter([mc_new],[b_value],marker='s',color='green')
# mlist = swarm_copy[ swarm_copy['magnitude']-mc_new>=0.0]['magnitude'] 
# (dmax,pval2nd) = scipy.stats.kstest( mlist, lambda x: 1.0-10**(-b_value*(x-mc_new)))
# ax4.plot([mc_new],[pval2nd],marker='o',color='green')


fig.tight_layout(w_pad=-0.5)
#fig.savefig('%s/b.png'%DIR_OUTPT_figs,dpi=75)
plt.show()

AttributeErrorTraceback (most recent call last)
<ipython-input-20-022e35ad1ed4> in <module>()
     30 
     31 #--- list of magnitudes
---> 32 mm_list = swarm_copy['magnitude'].to_list()
     33 mm_list = np.sort(mm_list)
     34 

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/pandas/core/generic.pyc in __getattr__(self, name)
   3612             if name in self._info_axis:
   3613                 return self[name]
-> 3614             return object.__getattribute__(self, name)
   3615 
   3616     def __setattr__(self, name, value):

AttributeError: 'Series' object has no attribute 'to_list'
In [21]:
##--- estimate b and mc
#
#
#with path.Path('/Users/Home/Desktop/Tmp/txt/myPYlibs'):
#    import plfit
#    import plvar
#    import plpva
#
##--- set min/max time to avoid temporal incompletenesss issue
#swarm_copy = swarm.copy()#DataFrameSubSet( swarm, 
#                          #   'date', 
#                           #  ( myTlohi.lo, 
#                            #   myTlohi.hi ) )
#
##--- power law fit clauset
#n_samp = 1000
#n_chunk = len( swarm_copy ) / n_samp
#
#i = 0
#fig = plt.figure(figsize=(8,4))
#ax = fig.add_subplot(111)
#ax.set_ylabel('m')
#ax2=ax.twinx()
#ax.scatter(swarm_copy['date'],swarm_copy['magnitude'],
#          alpha=0.1)
#
#fig2 = plt.figure(figsize=(8,4))
#ax21 = fig2.add_subplot(111)
#ax21.set_ylabel('m')
#ax21.scatter(swarm_copy['date'],swarm_copy['magnitude'],
#          alpha=0.1)
#
#MC = []
#B = []
#for j in xrange(n_chunk):
#    Mlist = 10**(swarm_copy.iloc[i:i+n_samp]['magnitude'])
#    fit = plfit.plfit( np.array(Mlist) ) #--- clauset
#    mc = log(fit[1],10) #--- completeness mag.
#    bval = fit[0]-1
#    MC.append(mc)
#    B.append(bval)
##    print mc, bval
#    ax2.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
#            [bval,bval],'-',color='red')
#
#    ax21.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
#            [mc,mc],'-',color='red')
#
#    i += n_samp
#Mlist = 10**(swarm_copy.iloc[i:i+len(swarm_copy)%n_samp]['magnitude'])
#fit = plfit.plfit( np.array(Mlist) ) #--- clauset
#mc = log(fit[1],10) #--- completeness mag.
#bval = fit[0]-1
#MC.append(mc)
#B.append(bval)
#
##---- plot bval
#ax2.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
#            [bval,bval],'-',color='red')
#ax2.tick_params(colors='red')
#ax2.set_ylabel('b',color='red')
#ax.set_xlabel('Date')
##print mc, bval
#fig.savefig('%s/b.png'%DIR_OUTPT_figs,dpi=75)
##--- plot mc
#ax21.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
#            [mc,mc],'-',color='red',label='mc')
#plt.legend()
#ax21.set_xlabel('Date')
#
##n, bins, patches = plt.hist(B, 50)
#
##--- error estimate
##[ error_b, error_mc, ntail ] = [ 0.0, 0.0, 0.0 ]
##[error_b, error_mc, ntail] = plvar.plvar(np.array(Mlist), 'silent')
##error_mc = log(1.0+error_mc/fit[1],10) #--- log(1+dM/M)
#plt.savefig('%s/mc.png'%DIR_OUTPT_figs,dpi=75)
In [22]:
# #--- goodness of the fit

# n_samp = 1000
# n_chunk = len( swarm_copy ) / n_samp
# ncols=4
# nrows = n_chunk / ncols + int(np.ceil(1.0*( len(swarm_copy) % n_samp) / n_samp))
# i = 0
# plt.figure(figsize=(ncols*4,nrows*4))
# for j in xrange(n_chunk):
#     Mlist = 10**(swarm_copy.iloc[i:i+n_samp]['magnitude'])
#     fit = plfit.plfit( np.array(Mlist) ) #--- clauset
#     mc = log(fit[1],10) #--- completeness mag.
#     bval = fit[0]-1
# #    [error_b, error_mc, ntail] = plvar.plvar(Mlist, 'silent') #--- error
#     plt.subplot(nrows,ncols,j+1)
#     Mlist=np.log10(Mlist)
#     nc = len(Mlist[Mlist>=mc])
#     error_b = (bval / nc / log(10.0) )**0.5
#     hist,bins = np.histogram(Mlist,bins=1024)
#     plt.plot(bins[:-1],-hist.cumsum()+len(Mlist))#,marker='o',markersize=6,linestyle='None')
#     xc=0.5*(bins[:-1]+bins[1:])
#     plt.plot(xc, len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r',label='$b=%2.1f\pm%3.2f$'%(bval,error_b))
#     plt.plot([mc,mc],[1,len(Mlist[Mlist>=mc])],label='$mc=%2.1f$'%(mc)) # len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r')
#     plt.yscale('log')
#     plt.xlim(swarm_copy['magnitude'].min(),swarm_copy['magnitude'].max())
#     plt.ylim(1,n_samp*10)
#     if j == 0:
#         plt.ylabel('N(size>m)')
#         plt.xlabel('m')
#     plt.legend()
#     i += n_samp

# Mlist = 10**(swarm_copy.iloc[i:i+len(swarm_copy)%n_samp]['magnitude'])
# fit = plfit.plfit( np.array(Mlist) ) #--- clauset
# mc = log(fit[1],10) #--- completeness mag.
# bval = fit[0]-1
# Mlist=np.log10(Mlist)
# nc = len(Mlist[Mlist>=mc])
# error_b = (bval / nc / log(10.0) )**0.5
# hist,bins = np.histogram(Mlist,bins=1024)

# plt.subplot(nrows,ncols,j+2)
# plt.plot(bins[:-1],-hist.cumsum()+len(Mlist))#,marker='o',markersize=6,linestyle='None')
# xc=0.5*(bins[:-1]+bins[1:])
# plt.plot(xc, len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r',label='$b=%2.1f\pm%3.2f$'%(bval,error_b))
# plt.plot([mc,mc],[1,len(Mlist[Mlist>=mc])],label='$m_c=%2.1f$'%(mc)) # len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r')
# plt.yscale('log')
# plt.xlim(swarm_copy['magnitude'].min(),swarm_copy['magnitude'].max())
# plt.ylim(1,n_samp)
# plt.legend()

# plt.savefig('%s/b_mult.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
    
In [23]:
# #--- plot complete catalog

# def ConvertDailyRate(hist, bin_edges ):
# #---convert data to daily rate     
#     t0 = pd.to_datetime( bin_edges[ 0 ] )
#     t1 = pd.to_datetime( bin_edges[ 1 ] )
#     delta_t = ( t1 - t0 ).total_seconds() / ( 60 * 60*24)
#     hist *= ( bin_edges[ 1 ] - bin_edges[ 0 ] ) / delta_t

# def ActivityRate( swarm ):
#     nbins = int( (swarm['date'].max()-swarm['date'].min()).days*4 ) #--- number of bins
    
#     tmax = swarm['date'].max().value #--- min/max
#     tmin = swarm['date'].min().value
    
#     hist, bin_edges = np.histogram(swarm['date'].apply(lambda x:x.value),                                   
#                     bins=np.linspace(tmin,tmax,nbins+1,endpoint=True),density=True) #--- probability dist.
#     hist *= len( swarm['date'] ) #--- int(hist).dt=n
#     cumm_number = np.cumsum(hist)*(bin_edges[1]-bin_edges[0]) #--- accumulated number
#     ConvertDailyRate( hist, bin_edges ) #--- daily activity
#     return bin_edges, hist, cumm_number
    
    
# #---------------------------------------------------------------------------------
# #-----------------
# #-----------------
# #-----------------
# #---------------------------------------------------------------------------------
   
# #--- completeness
# mc = 1.5



# #--- t0<t<t1
# #--- set min/max time to avoid temporal incompletenesss issue
# swarm_tmp = swarm.copy()#DataFrameSubSet( swarm, 
#                          #    'date', 
#                           #   ( swarm['date'].min(), 
#                            #    pd.to_datetime('2010-07-31') ) )
# #--- m > mc
# swarm_lohi = swarm_tmp.copy() #DataFrameSubSet( swarm_tmp, 
#                             # 'magnitude', 
#                              #( mc, sys.maxint ) ) 

# #--- spatial map
# #swarm_lohi.plot.scatter('longitude','latitude',
# #                        s=3**(swarm_lohi['magnitude']),
# #                        c='date',cmap='jet',
# #                        alpha = 0.4) #--- plot
    
# #--- temporal map
# #ax1.set_xlim(swarm_tmp['date'].min(),swarm_tmp['date'].max())
# for ii in xrange(10):
#     fig = plt.figure(figsize=(8,4))
#     ax1=fig.add_subplot(111)

#     #plt.xlabel('Time')
#     #plt.ylabel('M')

#     ax1.set_ylim(0,6)
#     tt=swarm[swarm['magnitude']==swarm['magnitude'].sort_values(ascending=False).iloc[ii]]['date'].iloc[0]
#     ax1.set_xlim(tt-datetime.timedelta(days=15),tt+datetime.timedelta(days=15))# pd.to_datetime('2014-12-04'),pd.to_datetime('2015-06-1'))
    
#     ax1.tick_params(axis='both',labelsize=18)
#     #ax1.tick_params(axis='x',rotation=18)
    
#     ax1.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
#                 s=2*np.exp(swarm_lohi['magnitude']),
#                 alpha=0.04,color='black')
        
#     #--- activity rate
#     bin_edges, hist, cumm_number = ActivityRate( swarm_lohi )
    
#     #--- plot
#     ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
#     ax2.plot(pd.to_datetime(bin_edges[:-1]),hist,'r-')
    
#     ax2.tick_params(axis='y',labelsize=18,colors='red')
    
#     ax2.set_yscale('log')
    
#     #ax2.set_ylim(1e1,1e3)
    
#     ax2.xaxis.set_major_formatter(dates.DateFormatter('%b %y'))
    
#     fig.autofmt_xdate()

# fig.savefig('%s/timeSeries_ok.png'%DIR_OUTPT_figs,bbox_inches='tight')
# plt.show()
In [24]:
#--- plot complete catalog

def ConvertDailyRate(hist, bin_edges ):
#---convert data to daily rate     
    t0 = pd.to_datetime( bin_edges[ 0 ] )
    t1 = pd.to_datetime( bin_edges[ 1 ] )
    delta_t = ( t1 - t0 ).total_seconds() / ( 60 * 60 )
    hist *= ( bin_edges[ 1 ] - bin_edges[ 0 ] ) / delta_t

def ActivityRate( swarm ):
    nbins = int( (swarm['date'].max()-swarm['date'].min()).days * 8 ) #--- number of bins
    
    tmax = swarm['date'].max().value #--- min/max
    tmin = swarm['date'].min().value
    
    hist, bin_edges = np.histogram(swarm['date'].apply(lambda x:x.value),                                   
                    bins=np.linspace(tmin,tmax,nbins+1,endpoint=True),density=True) #--- probability dist.
    hist *= len( swarm['date'] ) #--- int(hist).dt=n
    cumm_number = np.cumsum(hist)*(bin_edges[1]-bin_edges[0]) #--- accumulated number
    ConvertDailyRate( hist, bin_edges ) #--- daily activity
    return bin_edges, hist, cumm_number
    
    
#---------------------------------------------------------------------------------
#-----------------
#-----------------
#-----------------
#---------------------------------------------------------------------------------
   
#--- completeness
#mc = 1.5



#--- t0<t<t1
myTlohi = Tlohi( pd.to_datetime('2019-07-04'), 
                 pd.to_datetime('2019-07-16')) #swarm['date'].max())

#--- set min/max time to avoid temporal incompletenesss issue
swarm_tmp = DataFrameSubSet( swarm, 
                             'date', 
                             ( myTlohi.lo, 
                               myTlohi.hi ) )
#--- m > mc
swarm_lohi = swarm_tmp.copy() #DataFrameSubSet( swarm_tmp, 
                            # 'magnitude', 
                             #( mc, sys.maxint ) ) 

#--- spatial map
#swarm_lohi.plot.scatter('longitude','latitude',
#                        s=3**(swarm_lohi['magnitude']),
#                        c='date',cmap='jet',
#                        alpha = 0.4) #--- plot
    
#--- temporal map
fig = plt.figure(figsize=(8,3))#,dpi=150)
ax1=fig.add_subplot(111)

#plt.xlabel('Time')
#plt.ylabel('M')

ax1.set_ylim(0,7.2)
ax1.set_xlim(swarm_tmp['date'].min(),swarm_tmp['date'].max())

ax1.tick_params(axis='both',labelsize=18)
ax1.tick_params(axis='x',rotation=30)

ax1.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
            s=2*np.exp(swarm_lohi['magnitude']),
            alpha=0.04,color='black')

DrawFrame(ax1, (0.11,0.07),(0.32,0.1),0.01)

#--- activity rate
bin_edges, hist, cumm_number = ActivityRate( swarm_lohi )

#--- plot
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
ax2.plot(pd.to_datetime(bin_edges[:-1]),hist,'r-',linewidth=1)

ax2.tick_params(axis='y',labelsize=18,colors='red')

#ax2.set_yscale('log')

#ax2.set_ylim(1e1,3e2)

ax2.xaxis.set_major_formatter(dates.DateFormatter('%d %b %y'))

#fig.autofmt_xdate(bottom=.1)

#fig.savefig('%s/timeSeries_ok.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
plt.show()

TypeErrorTraceback (most recent call last)
<ipython-input-24-b697eb84d4a7> in <module>()
     68 ax1.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
     69             s=2*np.exp(swarm_lohi['magnitude']),
---> 70             alpha=0.04,color='black')
     71 
     72 DrawFrame(ax1, (0.11,0.07),(0.32,0.1),0.01)

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1715                     warnings.warn(msg % (label_namer, func.__name__),
   1716                                   RuntimeWarning, stacklevel=2)
-> 1717             return func(ax, *args, **kwargs)
   1718         pre_doc = inner.__doc__
   1719         if pre_doc is None:

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, **kwargs)
   4021             linewidths = rcParams['lines.linewidth']
   4022 
-> 4023         offsets = np.column_stack([x, y])
   4024 
   4025         collection = mcoll.PathCollection(

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/numpy/lib/shape_base.pyc in column_stack(tup)
    367             arr = array(arr, copy=False, subok=True, ndmin=2).T
    368         arrays.append(arr)
--> 369     return _nx.concatenate(arrays, 1)
    370 
    371 def dstack(tup):

TypeError: invalid type promotion
In [25]:
#--- evaluate fractal dimension

import geopy.distance
from math import *
import random as rnd

class bins:
    def __init__(self, nBin, xlo, xhi, ydim = 1, err = None):
        self.lo = xlo - 1.0e-10
        self.hi = xhi + 1.0e-10
        self.dx = (self.hi-self.lo)/nBin
        self.xlist = [0.0 for i in xrange(nBin)]
        self.kounter = [0 for i in xrange(nBin)]
        self.ylist = [[0.0 for j in xrange(ydim)] for i in xrange(nBin)]
        self.nBin = nBin
        self.err = err
        self.max_y = [[-sys.maxint for j in xrange(ydim)] for i in xrange(nBin) ]
        self.min_y = [[sys.maxint for j in xrange(ydim)] for i in xrange(nBin)]
        if err:
            self.ySQlist = [[0.0 for j in xrange(ydim)] for i in xrange(nBin)]
    def GetBin(self,x):
        return int(floor((x-self.lo)/self.dx))
    def whichBin(self,x,y, ibin=[] ):
        assert x >= self.lo, 'x=%s,self.lo=%s'%(10**x,10**self.lo)
        assert x < self.hi, 'x=%s,self.hi=%s'%(10**x,10**self.hi)
        nr = int(floor((x-self.lo)/self.dx))
        if ibin:
            ibin[ 0 ] = nr
        self.kounter[nr] += 1
        self.xlist[nr] += x
        for idim in xrange(len(y)):
            self.ylist[nr][idim] += y[idim]
            if y[idim] >= self.max_y[nr][idim]: #--- set max value
                self.max_y[nr][idim]=y[idim]
            if y[idim] <= self.min_y[nr][idim]:
                self.min_y[nr][idim]=y[idim]
            if self.err:
                self.ySQlist[nr][idim] += y[ idim ] * y[ idim ]

    def res(self, logScaleX = None, logScaleY = None, base = 10, SUM = None, MINMAX=None ):
         indices = xrange(10**6)
         someList = []
         for x,index in zip(self.xlist,indices):
             nb = self.kounter[index]
             if nb == 0: continue
             xbar = self.xlist[index]/nb
             ybar = [y/nb for y in self.ylist[index]]
             if self.err:
                 sigmaY = [ ( ysq / nb - YBAR * YBAR ) ** 0.5 / nb ** 0.5 for ysq, YBAR in zip( self.ySQlist[ index ], ybar )]
                 if SUM:
                     sigmaY = [ i * nb for i in sigmaY ]
             if SUM:
                 ybar = [y for y in self.ylist[index]]
             if MINMAX:
                 MAX_y = [y for y in self.max_y[index]]
             if logScaleX:
                 xbar = base ** xbar
             if logScaleY:
                 ybar = [ base ** item for item in ybar ]
             if self.err:
                 someList.append([ xbar, ybar, sigmaY ])
             elif MINMAX:
                 someList.append([ xbar, ybar, MAX_y ])
             else:
                 someList.append([ xbar, ybar ])
         return someList


class histogram( bins ):
    def res( self, Radial = None, logScale = None, normalize = True, base = 10.0, ACCUMLTD = None ):
        PDF = []
        self.nPoint = nPoint = sum( self.kounter )
        indices = xrange( sys.maxint )
        y_accm = nPoint
        for y, index in zip( self.kounter, indices ):
            if not ACCUMLTD and y == 0:
                continue
            if not y == 0:
                x = self.xlist[ index ] / y #self.lo + index * self.dx
            else:
                x = self.lo + index * self.dx
            Y = 1.0 * y
            dx = self.dx
            if logScale:
                x = base ** x
                dx = x * ( base ** self.dx - 1.0 )
#               print Y, dx
            if normalize:
                Y /= ( nPoint * dx )
                if Radial:
                    Y /= ( 2.0 * pi * x )
#           PDF.append( [ x, Y ] )
#           PDF.append( [ x + dx, Y ] )
#           PDF.append( [ x + 0.5 * dx, Y, 0.0, ( 1.0 * y_accm / nPoint if normalize else 1.0 * y_accm )  ] )
            error_std = 0.0
            if self.err:
                error_std = sqrt( nPoint * Y * dx ) / ( nPoint * dx ) #--- poisson
                error_std = sqrt( nPoint * Y * dx * ( 1.0 - Y * dx ) ) / ( nPoint * dx ) #--- bi-nomial
            PDF.append( [ x, Y, 0.0, ( 1.0 * y_accm / nPoint if normalize else 1.0 * y_accm ), error_std ] )
            y_accm -= y
        return PDF

def GetCartesian( dff ):
    df = dff.copy()
    xlo = df['longitude'].min()
    xhi = df['longitude'].max()
    ylo = df['latitude'].min()
    yhi = df['latitude'].max()
    getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
    getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
    df[ 'x(km)' ] = df[ 'longitude' ].apply( getDistX ) 
    df[ 'y(km)' ] = df[ 'latitude' ].apply( getDistY ) 
    df[ 'z(km)' ] = df['depth']
    return df

def fractalDimension2nd( coord ):
    #--- sort
    if type( coord ) == type( [] ):
        coord = ListToDict( coord )
    points = coord.keys()
    points.sort()
    hsObject = histogram( 18 * 8, log( 1e-10,10), log( 1e8, 10 ) )
    for point_i in points:
        for point_j in points:
            if not point_i < point_j: #--- pairs ij with i<j
                continue
            rij = sum( [ ( i - j ) ** 2 for i, j in zip( coord[ point_i ], coord[ point_j ] ) ] ) # ** 0.5
            assert rij > 0, 'rij=%s,coord[ %s ]=%s, coord[ %s ]=%s' %(rij,point_i,coord[ point_i ], point_j, coord[ point_j ] )
            hsObject.whichBin( 0.5 * log( rij,10 ), [ 1.0 ] )
    for items in hsObject.res( logScale = True, normalize = True, ACCUMLTD = True ):
        if items[ 3 ] > 0.0: 
            yield items[ 0 ], items[ 3 ]

#---------------------------------------------------------------------------------
#-----------------
#-----------------
#-----------------
#---------------------------------------------------------------------------------
mc = 2.0

#--------------------
#----- subset
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi = swarm_lohi.sample( n = 2000 ) #--- sample

#--------------------
#--- cartesian coords
#--------------------
swarm_cartesian = GetCartesian( swarm_lohi )

#--------------------
#--- evaluate df
#--------------------
tmp_coord = swarm_cartesian['x(km)'].apply(lambda x: [x+rnd.random()*1e-6]) +\
            swarm_cartesian['y(km)'].apply(lambda x: [x+rnd.random()*1e-6]) +\
            swarm_cartesian['z(km)'].apply(lambda x: [x+rnd.random()*1e-6])
tmp_coord = tmp_coord.to_dict()
dict_NR = fractalDimension2nd( tmp_coord ) #, dmin = 1.0e-02 )

#--------------------
#--- scattered plot
#--------------------
swarm_cartesian.plot.scatter('x(km)','y(km)',
                        s=3**(swarm_lohi['magnitude']),
                        c='date',cmap='jet',
                        alpha = 0.4) 
plt.savefig('%s/map2d.png'%DIR_OUTPT_figs,bbox_inches='tight')
#--------------------
#--- N(r) vs r
#--------------------
plt.figure( figsize = (4,4))
plt.xlabel('r(km)')
plt.ylabel('N(r)')
plt.xlim(1e-3,1e2)
plt.ylim(1e-8,1)
plt.yscale('log')
plt.xscale('log')
d_f = 2.2 #1.7
for i in dict_NR:
    plt.plot([i[ 0 ]],
             [1-i[ 1 ]],
             'o',color='black') #--- accumulated
    plt.plot(i[ 0 ],
             (1-i[ 1 ])/i[0]**d_f,
             '.',color='red') #--- accumulated
plt.show()
In [26]:
#--- evaluate fractal dimension (correltion dimension )

def getRmat2d( df_complete ):
    nmax = len( df_complete )
    rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    for i in xrange( nmax ):
        #--- distance matrix
        df_dx = df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ]
        df_dy = df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ]
        df_sq = ( df_dx*df_dx+df_dy*df_dy ) ** 0.5
        df_sq[ : i ] = 0
        rmat[ i ] = np.array(df_sq)
    return np.array( rmat ) 

def getRmat3d( df_complete ):
    nmax = len( df_complete )
    rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    for i in xrange( nmax ):
        #--- distance matrix
        df_dx = df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ]
        df_dy = df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ]
        df_dz = df_complete[ 'z(km)' ] - df_complete[ 'z(km)' ][ i ]
        df_sq = ( df_dx*df_dx+df_dy*df_dy+df_dz*df_dz ) ** 0.5
        df_sq[ : i ] = 0
        rmat[ i ] = np.array(df_sq)
    return np.array( rmat ) 

#def GetCartesian( dff ):
#    df = dff.copy()
#    xlo = df['longitude'].min()
#    xhi = df['longitude'].max()
#    ylo = df['latitude'].min()
#    yhi = df['latitude'].max()
#    getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
#    getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
#    df[ 'r(km)' ] = df[ 'longitude' ].apply( getDistX ) + df[ 'latitude' ].apply( getDistY ) * 1j
#    return df


# --------------------
# ----- subset
# --------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 

#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )

#swarm_lohi = DataFrameSubSet( swarm_lohi, 
 #                            'x(km)', 
  #                           ( 180,300)) 
#swarm_lohi = DataFrameSubSet( swarm_lohi, 
 #                            'y(km)', 
  #                           ( 100,250)) 

                              
#--- we need the distance matrix! 
#rmat2d = getRmat2d( swarm_lohi )
rmat3d = getRmat3d( swarm_lohi )
#r_list2d = rmat2d[rmat2d.nonzero()]
r_list3d = rmat3d[rmat3d.nonzero()]

# #--- histogram
r_min = log(np.min(r_list3d),10) #min(log(np.min(r_list2d),10),log(np.min(r_list3d),10))
r_max = log(np.max(r_list3d),10) #max(log(np.max(r_list2d),10),log(np.max(r_list3d),10))
nbins = int(ceil(r_max - r_min))*4
#2d
hist3d, edges = np.histogram(r_list3d,
                         bins=np.logspace(r_min, r_max, nbins) )
r_bins, edges = np.histogram(r_list3d,
                         bins=np.logspace(r_min, r_max, nbins), weights = r_list3d )

r_bins3d = r_bins[hist3d != 0]
hist3d = hist3d[hist3d != 0]
r_bins3d /= hist3d

# #--- normalize
hist_acc_3d = np.cumsum(hist3d)*1.0
hist_acc_3d /= hist_acc_3d[-1]

# #3d
# hist3d, edges = np.histogram(r_list3d,
#                          bins=np.logspace(r_min, r_max, nbins) )
# r_bins, edges = np.histogram(r_list3d,
#                          bins=np.logspace(r_min, r_max, nbins), weights = r_list3d )

# r_bins3d = r_bins[hist3d != 0]
# hist3d = hist3d[hist3d != 0]
# r_bins3d /= hist3d

# #--- normalize
# hist_acc_3d = np.cumsum(hist3d)*1.0
# hist_acc_3d /= hist_acc_3d[-1]

##--------------------
##--- scattered plot
##--------------------
#plt.scatter(swarm_lohi['x(km)'],
#            swarm_lohi['y(km)'],
#                        s=3**(swarm_lohi['magnitude']),
#                        cmap='jet',
#                        alpha = 0.1) 
#                             
##--- plot
##d_f = 1.4
#fig = plt.figure(figsize=(4,4))
#ax = fig.add_subplot(111)
#ax.set_xlim(10**r_min,10**r_max)
#ax.set_xlabel('$r(km)$')
#ax.set_ylabel('$N(r)$')
#ax.set_xscale('log')
#ax.set_yscale('log')
#ax.errorbar( r_bins, hist_acc, 
#            yerr = hist_acc/hist**.5,
#            color='black',marker='o')
#
#ax.errorbar( r_bins, 10000*hist_acc/r_bins**d_f,  #--- rescaled
#            yerr = hist_acc/hist**.5/r_bins**d_f,
#            color='red',marker='o',label='$N/r^{%s}$'%d_f)
#ax.legend()
#
#fig.savefig('%s/corrDim.png'%DIR_OUTPT_figs,bbox_inches='tight')
In [27]:
#--- plot
d_f = 2.1 #1.7
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.axis([1e-5,1e3,1e-8,1e0])
ax.set_xlabel('$r(km)$')
ax.set_ylabel('$N(r)$')
ax.set_xscale('log')
ax.set_yscale('log')
ax.errorbar( r_bins3d, hist_acc_3d, 
            yerr = hist_acc_3d/hist3d**.5,
            color='black',marker='o')
# ax.errorbar( r_bins3d, hist_acc_3d, 
#             yerr = hist_acc_3d/hist3d**.5,
#             color='red',marker='s')

#ax.errorbar( r_bins, 10000*hist_acc/r_bins**d_f,  #--- rescaled
#            yerr = hist_acc/hist**.5/r_bins**d_f,
#            color='red',marker='o',label='$N/r^{%s}$'%d_f)
ax.plot(r_bins3d, r_bins3d**d_f,
        '-.r',label='$r^{%s}$'%d_f)

ax.legend()

fig.savefig('%s/corrDim.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()
In [28]:
#--- evaluate fractal dimension (box counting )

import geopy.distance
from math import *
import random as rnd
#import warnings

#warnings.filterwarnings('ignore') #--- get rid of warnings

#def GetCartesian( dff ):
#    df = dff.copy()
#    xlo = df['longitude'].min()
#    xhi = df['longitude'].max()
#    ylo = df['latitude'].min()
#    yhi = df['latitude'].max()
#    getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
#    getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
#    df[ 'r(km)' ] = df[ 'longitude' ].apply( getDistX ) + df[ 'latitude' ].apply( getDistY ) * 1j
#    df[ 'x(km)' ] = df[ 'longitude' ].apply( getDistX )
#    df[ 'y(km)' ] = df[ 'latitude' ].apply( getDistY )
#    return df

#--------------------
#----- subset
#--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# #swarm_lohi = swarm_lohi.sample( n = 1000 ) #--- sample


# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )
# #swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
# #                             'x(km)', 
# #                             ( 200, 1000 ) ) 
# #swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
# #                             'y(km)', 
# #                             ( -1000, 220 ) ) 

# #--------------------
# #--- scattered plot
# #--------------------
# plt.figure(figsize=(4,4))
# #plt.xlim(0,np.max(r_list))
# #plt.ylim(-np.max(r_list)/2,np.max(r_list)/2)
# plt.xlabel('x(km)')
# plt.ylabel('z(km)')
# plt.scatter(swarm_lohi['x(km)'],
#             swarm_lohi['z(km)'],
#                         s=3**(swarm_lohi['magnitude']),
#                         cmap='jet',
#                         alpha = 0.01) 
# #--------------------
# #--- box counting
# #--------------------
# #tmp=pd.DataFrame({'x(m)':swarm_lohi['y(m)'],'y(m)':swarm_lohi['z(m)']}) #,'z(m)':swarm_lohi['z(m)']})
# #rmat = getRmat2d( tmp )
# #rmat = getRmat2d( swarm_lohi )

# #r_list = rmat[rmat.nonzero()]
# #r_min = np.min(r_list)
# #r_max = np.max(r_list)
# #
# plt.savefig('%s/map2d.png'%DIR_OUTPT_figs,bbox_inches='tight')
In [29]:
# #--- evaluate fractal dimension (box counting )


# def fractalDimension3d( x, y, z, L, dmin = None ):
#     x = np.array(x)
#     y=np.array(y)
#     z=np.array(z)
# #    xlo = x.min()
# #    xhi = x.max()
# #    ylo = y.min()
# #    yhi = y.max()
# #    zlo = z.min()
# #    zhi = z.max()
# #    L = max( zhi - zlo, yhi - ylo, xhi - xlo )
#     #---
#     xc = x.mean() #--- center
#     yc = y.mean()
#     zc = z.mean()
    
#     xo = xc - L #--- origin
#     yo = yc - L
#     zo = zc - L
# #    L *= 2.0
#     #---
#     kounter = 0  #--- number of divisions
#     dx = dy = dz = L  #--- scale
#     dl = []
#     N = []
#     while dx - dmin > 0.0: #--- do until dx < dmin
#         d = {}
#         index_i = np.array( map(int,np.floor( ( y - yo ) / dy ) ) )
#         index_j = np.array( map( int, np.floor( ( x - xo ) / dx ) ) )
#         index_k = np.array( map( int, np.floor( ( z - zo ) / dz ) ) )
        
#         inside_cell_lo = np.all([index_i>=0,index_j>=0,index_k>=0.0],axis=0) #--- point is inside?
#         inside_cell_hi = np.all([index_i<2**kounter,index_j<2**kounter,index_k<2**kounter],axis=0) #
#         inside_cell = np.all([inside_cell_lo,inside_cell_hi],axis = 0)

        
#         cell_id = index_k * 4 ** kounter + index_i * ( 2**kounter ) + index_j
#         #---
#         dl+=[dx]
#         N+=[len( set( cell_id ) )]
#         #---
#         kounter += 1
#         dx = dy = dz = L / 2**kounter
#     assert dx - dmin <= 0.0
#     #---
#     return np.array(dl), np.array(N)


# def fractalDimension2d( x, y, L, dmin = None ):
#     x = np.array(x)
#     y=np.array(y)
# #    xlo = x.min()
# #    xhi = x.max()
# #    ylo = y.min()
# #    yhi = y.max()
# #    L = max( yhi - ylo, xhi - xlo )
#     #---
#     xc = x.mean() #--- center
#     yc = y.mean()

#     xo = xc - L #--- origin
#     yo = yc - L
# #    L *= 2.0
#     #---
#     kounter = 0  #--- number of divisions
#     dx = dy = L / 2**kounter  #--- scale
#     dl = []
#     N = []
# #    pdb.set_trace()
# #    fig = plt.figure(figsize=(4,4))
# #    ax=fig.add_subplot(111)
# #    plt.scatter(x,y)
# #    ax.add_patch( patches.Rectangle(xy=(xo,yo), width=L, 
# #                                    height=L, linewidth=1,
# #                                    clip_on=False,facecolor=None,edgecolor='black',fill=None) ) 
    
#     while dx - dmin > 0.0: #--- do until dx < dmin
#         d = {}
        
#         index_i = np.array( map(int,np.floor( ( y - yo ) / dy ) ) ) #--- assingn index
#         index_j = np.array( map( int, np.floor( ( x - xo ) / dx ) ) )

#         inside_cell_lo = np.all([index_i>=0,index_j>=0],axis=0) #--- point is inside?
#         inside_cell_hi = np.all([index_i<2**kounter,index_j<2**kounter],axis=0) #
#         inside_cell = np.all([inside_cell_lo,inside_cell_hi],axis = 0)
        
#         index_i = index_i[ inside_cell ] #--- subset of the index array
#         index_j = index_j[ inside_cell ]

#         cell_id = index_i * ( 2**kounter ) + index_j
#         #---
#         dl+=[dx]
#         N+=[len( set( cell_id ) )]
#         #---
#         kounter += 1
#         dx = dy = L / 2**kounter
#     assert dx - dmin <= 0.0
#     #---
#     return np.array(dl), np.array(N)


# fig = plt.figure(figsize=(10,4))
# ax = fig.add_subplot(121)
# ax2 = fig.add_subplot(122)
# ax.set_xlabel('$r(km)$')
# ax.set_ylabel('$N(r)$')
# ax.set_xscale('log')
# ax.set_yscale('log')
# ax2.set_xlabel('$r(km)$')
# ax2.set_ylabel(r'$N(r)\times r^{%s}$'%d_f)
# ax2.set_xscale('log')
# ax2.set_yscale('log')
# dl={}
# N={}
# count = 0
# r_min = np.min(r_list)
# r_max = np.max(r_list)
# for L in np.logspace(log(r_min,2),log(r_max,2),num=6,base=2.0):
# #for L in [2*420.0]:
# #    dl[count], N[count] = fractalDimension2d( swarm_lohi['x(km)'],
# #                                              swarm_lohi['y(km)'], L, dmin=r_min ) # , swarm_lohi['z(m)'], dmin=r_min )
#     dl[count], N[count] = fractalDimension3d( swarm_lohi['x(km)'],
#                                              swarm_lohi['y(km)'], 
#                                              swarm_lohi['z(km)'], 
#                                              L, dmin=r_min ) 

# ##--- try random
# #x=np.array([rnd.random() for i in xrange(1000)])
# #y=0.2*x #[rnd.random() for i in xrange(10000)]
# ##z=x #[rnd.random() for i in xrange(10000)]
# #tmp=pd.DataFrame({'x(m)':x,'y(m)':y})#,'z(m)':z})
# #rmat = getRmat2d( tmp )
# #r_list = rmat[rmat.nonzero()]
# #r_min = min(r_list)
# #sdict = fractalDimension2d( x,y,
# #                           dmin=r_min )
# #plt.scatter(x,y)

# #--- plot
#     #ax.set_xlim(10**r_min,10**r_max)
#     #for key in sdict:
#     ax.plot(dl[count],N[count],'-o')
#     ax2.plot(dl[count],N[count]*dl[count]**d_f,marker='o',label='$L=%2.1e$'%L)#,linestyle='None')
#     count+=1
# ax2.legend(bbox_to_anchor=(1, 1))
# fig.tight_layout()

# fig.savefig('%s/boxCount.png'%DIR_OUTPT_figs,bbox_inches='tight')
    
In [30]:
Df = 2.1 #1.7 #2.2 #2.0 #1.7
bval = 0.8
mc=2.0 #1.5
In [31]:
#--- trig analysis (new) (it's a large data sets. run the simulation on a cluster and use the outputs)

from scipy.sparse import lil_matrix
import time
from IPython.display import display
import datetime
import pdb
import numpy.linalg
import numpy as np

def UpperTriang(a):
    il1 = np.tril_indices(a.shape[0])
    a[il1] = 0
    return a
    
def getTmat( df_complete ):
    nmax = len( df_complete )
    prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
    tmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    tmat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    for i in xrange( nmax ):
        df_dt = ( df_complete[ 'date' ] - df_complete[ 'date' ][ i ] ).dt.total_seconds() * prefact #--- time diff between i-th event and all subsequent events	
        tmat[ i ] = np.array(df_dt)
        df_dt[ : i ] = df_dt[ i ] #--- must have an upper triangular matrix 
        tmat_ut[ i ] = np.array(df_dt)

        #---
    return np.array( tmat ), np.array( tmat_ut )


def AddPositionList(df_complete):
    df_complete['r(m)'] = df_complete['x(m)'].apply(lambda x: [x])+ \
                          df_complete['y(m)'].apply(lambda x: [x])+ \
                          df_complete['z(m)'].apply(lambda x: [x])

def getRmat( df_complete ):
    nmax = len( df_complete )
    rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    rmat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    for i in xrange( nmax ):
        #--- distance matrix
        df_dx = ( df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ] )
        df_dy = ( df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ] )
        df_dz = ( df_complete[ 'z(km)' ] - df_complete[ 'z(km)' ][ i ] ) #comment if 2-d
        df_sq = ( df_dx * df_dx + df_dy * df_dy + df_dz * df_dz ) ** 0.5
        rmat[ i ] = np.array(df_sq)
        df_sq[ : i ] = 0
        rmat_ut[ i ] = np.array(df_sq)
    return np.array( rmat ), np.array( rmat_ut )

def getMmat( df_complete ):
    nmax = len( df_complete )
    m_mat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    m_mat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
    for i in xrange( nmax ):
        #--- magnitude
        df_m = pd.Series( [ df_complete[ 'magnitude' ][ i ] ] * nmax )
        m_mat[ i ] = np.array(df_m)
        df_m[ : i ] = 0
        m_mat_ut[ i ] = np.array(df_m)
    return np.array( m_mat ), np.array( m_mat_ut ) 

def vectorizedAnalysis( df_complete ):
    #--- setup t, r, m
    t0 = time.time()
    tmat0, tmat = getTmat( df_complete )
    print 'setting up tmat:%s s'%(time.time() - t0)
    t0 = time.time()
    r_mat0, r_mat = getRmat( df_complete )
    print 'setting up rmat:%s s'%(time.time() - t0)
    t0 = time.time()
    m_mat0, m_mat = getMmat( df_complete )
    print 'setting up m_mat:%s s'%(time.time() - t0)
      
        
#    #--- make upper-triangular matrices    
#    t0 = time.time()
#    tmat = UpperTriang(np.copy(tmat0))
#    r_mat = UpperTriang(np.copy(r_mat0))
#    m_mat = UpperTriang(np.copy(m_mat0))
#    print 'setting up ut matrix:%s s'%(time.time() - t0)
    
    #--- nij
    t0 = time.time()
    NIJ = tmat * r_mat ** ( Df ) * 10 ** ( - bval * ( m_mat - mc ) )
    TIJ = tmat * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
    RIJ = NIJ / TIJ #r_mat ** ( Df ) * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
    print 'setting up NIJ:%s s'%(time.time() - t0)
    nmax = len( df_complete )
    N_sparse = lil_matrix( ( nmax, nmax ) ) #--- sparse matrix
    T_sparse = lil_matrix( ( nmax, nmax ) )
    R_sparse = lil_matrix( ( nmax, nmax ) )
    for junk, j in zip( NIJ, xrange( sys.maxint ) ): 
        if j == 0:
            continue
        x = min( NIJ[ : j, j ][ NIJ[ : j, j ] != 0 ] ) #--- min nij
        assert x != 0.0
        rowi = np.where( NIJ[ : j, j ] == x )[ 0 ][ 0 ] #--- find row
        assert rowi <= j
        N_sparse[ rowi, j ] = x #--- insert
        T_sparse[ rowi, j ] = TIJ[ rowi, j ]
        R_sparse[ rowi, j ] = RIJ[ rowi, j ]    
    return N_sparse, (T_sparse, tmat0 ), ( R_sparse, r_mat0 ), m_mat0


#--------------------
#----- subset
#--------------------
#--- t0<t<t1
# myTlohi = Tlohi( pd.to_datetime('2019-07-04'), 
#                  swarm['date'].max())

# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'date', 
#                              ( myTlohi.lo, myTlohi.hi ) ) 

swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)

#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )

#AddPositionList( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)

#--- vectorized
t0 = time.time()
N_sparse, (T_sparse,tmat), (R_sparse,rmat), m_mat = vectorizedAnalysis( swarm_lohi )
print 'duration=%s s'%(time.time() - t0)
setting up tmat:4.89720392227 s
setting up rmat:3.67028999329 s
setting up m_mat:3.49758481979 s
setting up NIJ:3.08642578125 s
duration=15.8014199734 s
In [32]:
#--- random analysis

def shuffle(amat):
    ind_list = range(amat.shape[0])
    np.random.shuffle(ind_list)
    amat_row = amat[ind_list]
    return amat_row[:,ind_list]

def GetNij(rmat,tmat,m_mat,df_complete):
    
    tmat = shuffle(np.abs(tmat))
    r_mat = shuffle(rmat)
    m_mat = shuffle(m_mat)
    
    #--- nij
    NIJ = tmat * r_mat ** ( Df ) * 10 ** ( - bval * ( m_mat - mc ) )
    TIJ = tmat * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
    RIJ = r_mat ** ( Df ) * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
    nmax = len( df_complete )
    N_sparse = lil_matrix( ( nmax, nmax ) ) #--- sparse matrix
    T_sparse = lil_matrix( ( nmax, nmax ) )
    R_sparse = lil_matrix( ( nmax, nmax ) )
    for junk, j in zip( NIJ, xrange( sys.maxint ) ): 
        if j == 0:
            continue
            
        x = min( NIJ[ : j, j ][ NIJ[ : j, j ] != 0 ] ) #--- min nij
        assert x != 0.0
        rowi = np.where( NIJ[ : j, j ] == x )[ 0 ][ 0 ] #--- find row
        N_sparse[ rowi, j ] = x #--- insert
        T_sparse[ rowi, j ] = TIJ[ rowi, j ]
        R_sparse[ rowi, j ] = RIJ[ rowi, j ]    
    return N_sparse, (T_sparse, tmat ), ( R_sparse, r_mat), m_mat

N_sparse_rnd = {}
T_sparse_rnd = {}
R_sparse_rnd = {}
tmat_rnd = {}
rmat_rnd = {}
nsamp = 1
tmat = np.abs(tmat)
for i in xrange(nsamp):
    print i
    N_sparse_rnd[i], (T_sparse_rnd[i],tmat_rnd[i]), (R_sparse_rnd[i],rmat_rnd[i]), junk = GetNij(rmat,tmat,m_mat,swarm_lohi)
0
In [33]:
# #--- save sparse matrices

# import scipy.sparse

# for mat, title in zip([N_sparse,T_sparse,R_sparse],
#                      ['N_sparse_matrix.npz','T_sparse_matrix.npz','R_sparse_matrix']):
#     scipy.sparse.save_npz('%s/%s'%(DIR_OUTPT,title), mat.tocsr()) #--- output

# for mat, title in zip([N_sparse_rnd[0],T_sparse_rnd[0],R_sparse_rnd[0]],
#                      ['N_sparse_rnd_matrix.npz','T_sparse_rnd_matrix.npz','R_sparse_rnd_matrix']):
#     scipy.sparse.save_npz('%s/%s'%(DIR_OUTPT,title), mat.tocsr()) #--- output

# #swarm_shuffled.to_csv('%s/swarm_shuffled.csv'%DIR_OUTPT)
In [34]:
# #--- load sparse matrices
# #--- 1- run on the cluster: 
# #--- 2- sbatch --mem=8gb --partition=single  --time=02:59:59 -n 1 ./oarScript.sh
# #--- 3- oarScript.sh: source activate conda-env; python ./swarmEMC.py
# #--- 4- copy: scp arc:/home/kamran.karimi1/Project/seismic/dataset/El_Mayor_Cucpah/* ./swarm/dataset/El_Mayor_Cucpah/

# import scipy.sparse
# DIR_OUTPT0 = './runRidgeCrest/Run0'
# N_sparse = scipy.sparse.load_npz('%s/N_sparse_matrix.npz'%DIR_OUTPT0 )
# T_sparse = scipy.sparse.load_npz('%s/T_sparse_matrix.npz'%DIR_OUTPT0 )
# R_sparse = scipy.sparse.load_npz('%s/R_sparse_matrix.npz'%DIR_OUTPT0 )
# #rmat = np.load('%s/rmat.npy'%DIR_OUTPT)
# #tmat = np.load('%s/tmat.npy'%DIR_OUTPT)

# N_sparse_rnd = {}
# T_sparse_rnd = {}
# R_sparse_rnd = {}

# nsamp = 1
# N_sparse_rnd[0] = scipy.sparse.load_npz('%s/N_sparse_rnd_matrix.npz'%DIR_OUTPT0 )
# T_sparse_rnd[0] = scipy.sparse.load_npz('%s/T_sparse_rnd_matrix.npz'%DIR_OUTPT0 )
# R_sparse_rnd[0] = scipy.sparse.load_npz('%s/R_sparse_rnd_matrix.npz'%DIR_OUTPT0 )
# #rmat_rnd = np.load('%s/rmat_rnd.npy'%DIR_OUTPT)
# #tmat_rnd = np.load('%s/tmat_rnd.npy'%DIR_OUTPT)

# #swarm_shuffled = pd.read_csv('%s/swarm_shuffled.csv'%DIR_OUTPT)
# #swarm_shuffled_incomplete = pd.read_csv('%s/swarm_shuffled_incomplete.csv'%DIR_OUTPT)
In [35]:
##--- matrix map of nij & nij_rand
#
#plt.subplot(2,1,1)
#plt.pcolormesh(np.log(N_sparse.toarray()),cmap='jet')
#plt.subplot(2,1,2)
#plt.pcolormesh(np.log(N_sparse_rnd.toarray()),cmap='jet')
In [36]:
# mpl.rcParams.update(mpl.rcParamsDefault)
# warnings.filterwarnings('ignore') #--- get rid of warnings

# rc('text', usetex=True)
# font = {'size'   : 26}
# matplotlib.rc('font', **font)

#---------------------------
#--- compute threshold n_th
#---------------------------

def PLOT(edges, hist, cdf, title ):
    #--- plot
    plt.figure(figsize=(8,4))
    plt.subplot(1,2,1)
#    plt.xlim(1e-13,1)
#    plt.ylim(1e-1,1e13)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('nij')
    plt.ylabel('P(nij)')
    plt.title( title )
    plt.plot(edges[:-1],hist,'.-')
    
    plt.subplot(1,2,2)
#    plt.xlim(1e-13,1)
#    plt.ylim(1e-3,1)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('nij')
    plt.ylabel('cdf')
    plt.plot(edges[:-1],cdf,'.-')
    plt.plot([n_thresh for i in edges[:-1]], cdf)
#---    
def HIST( nij, nbins ):
    #--- histograms
    nij_TxR = nij['T*']*nij['R*']
    hist, edges = np.histogram( nij_TxR,
                            bins=np.logspace(np.log(nij_TxR.min()),np.log(nij_TxR.max()),nbins),
                            density=True)
    #--- accumulated
    cdf, edges = np.histogram( nij_TxR,
                           bins=np.logspace(np.log(nij_TxR.min()),np.log(nij_TxR.max()),nbins,endpoint=True))
    cdf=np.cumsum(cdf)
    #assert cdf[-1] == len(nij_TxR), '%s,%s'%(cdf[-1], len(nij_TxR))
    N = 1.0 * len(nij_TxR)
    cdf = cdf / N
    return edges, hist, cdf

def reset_xticks(nn,ax,(TMIN,TMAX)):
    ax.set_xticks( np.linspace(ax.axis()[0],ax.axis()[1],nn) ) 
    labels=np.linspace(TMIN,TMAX,nn)
    labels = [r'$10^{%i}$'%i for i in labels]
    ax.set_xticklabels(labels)

def reset_yticks(nn,ax,(RMIN,RMAX)):
    ax.set_yticks( np.linspace(ax.axis()[2],ax.axis()[3],nn) ) 
    labels=np.linspace(RMIN,RMAX,nn)
    labels = [r'$10^{%i}$'%i for i in labels]
    ax.set_yticklabels(labels)
                 

def GetTrigListLoHi( nij_trig, col, (tlo,thi), swarm_lohi ):
    list_of_mothers = nij_trig.groupby(by='Parent_id').groups.keys() 
    list_of_mothers.sort()
    tmp_df = swarm_lohi.loc[ list_of_mothers ]
    tmp_df = tmp_df[ (tmp_df[ col ] >= tlo) & 
                 (tmp_df[ col ] < thi) ]
    return nij_trig[ pd.DataFrame([nij_trig['Parent_id'] == i for i in tmp_df.index]).any() ]

def GetInterpolatedData( df0, (TMIN,TMAX),(RMIN,RMAX) ):
    nbins_per_decade = 16
    nbins_x=nbins_y=int(TMAX-TMIN)*nbins_per_decade
    df=df0.copy()

    df['T*']=df[ 'T*' ].apply(lambda x: log(x,10))
    df['R*']=df[ 'R*' ].apply(lambda x: log(x,10))
    heatmap, xedges, yedges = np.histogram2d( df[ 'T*' ], df[ 'R*'], bins=[np.linspace(TMIN,TMAX,nbins_x+1), np.linspace(RMIN,RMAX,nbins_y+1)], normed=True)
    heatmap *= len( df )
    heatmap = gaussian_filter( heatmap, sigma = nbins_per_decade/4 )
    return heatmap

def GetDistance( df1, df2 ):
    df_sq = (np.array( df1['x(km)'] ) - np.array( df2['x(km)'] ) ) ** 2 + \
            (np.array( df1['y(km)'] ) - np.array( df2['y(km)'] ) ) ** 2 + \
            (np.array( df1['z(km)'] ) - np.array( df2['z(km)'] ) ) ** 2 
    return pd.Series( df_sq ** 0.5 )

def RemovePair( dff, cs, catalog ):
    pref = 1.0 / 3600.0 / 24 #--- daily
    df = dff.copy()
    
    series_r = GetDistance(catalog.loc[df['Parent_id']],catalog.loc[df['Event_id']])
     
#    series_r = pd.Series(np.array( catalog.loc[df['Parent_id']]['r(km)'] ) -\
#                         np.array( catalog.loc[df['Event_id']]['r(km)'] ) )
    assert len(series_r[series_r==np.nan]) == 0
    df['R'] = series_r.copy() #abs()

    series_t = pd.Series(np.array( catalog.loc[df['Event_id']]['date'] ) -\
                         np.array( catalog.loc[df['Parent_id']]['date'] ) )
    df['T'] = series_t.apply(lambda x:x.total_seconds()*pref)
    assert len ( df[ df['R'] == 0.0 ] ) == 0, '%s'%display( df[ df['R'] == 0.0 ] )
    assert len ( df[ df['T'] == 0.0 ] ) == 0, '%s'%display( df[ df['T'] == 0.0 ] )
    return df [ df[ 'R' ] <= df[ 'T' ] * cs ], df [ df[ 'R' ] > df[ 'T' ] * cs ]


#--------------------
#----- set parameters
#--------------------
#--- set quantile
#quantile = 0.5 #0.35 #--- comment if passed by argument
#--- wave speed
cs = 3 * (24*3600) #--- km per day
#--- cut-off time
dt_cut = 0.0 #--- day
#--- plot scattered
# TMIN=-5
# TMAX=3
# RMIN=-4
# RMAX=4


#--------------------
#----- subset (the same as trig. analysis)
#--------------------
#--- t0<t<t1
# myTlohi = Tlohi( pd.to_datetime('2019-07-04'), 
#                  swarm['date'].max())

# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'date', 
#                              ( myTlohi.lo, myTlohi.hi ) ) 

swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)

#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )

#AddPositionList( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)


#--------------------
#--- store actual data in a dataframe
#--------------------
rows, cols = N_sparse.nonzero()
tlist = T_sparse[rows, cols].toarray()[0]
rlist = R_sparse[rows, cols].toarray()[0]
nij = pd.DataFrame({'Event_id':cols, 'Parent_id':rows,'T*':tlist,'R*':rlist})

#--------------------
#--- store random data in a dataframe
#--------------------
nij_rnd = {}
for i in xrange(nsamp):
    rows, cols = N_sparse_rnd[i].nonzero()
    tlist = T_sparse_rnd[i][rows, cols].toarray()[0]
    rlist = R_sparse_rnd[i][rows, cols].toarray()[0]
    nnn = len(cols)
    nij_rnd[i] = pd.DataFrame({'Event_id':cols, 'Parent_id':rows,'T*':tlist,'R*':rlist},index=range(i*nnn,(i+1)*nnn))


#--------------------------------------------
#----- only include mother events with 
#----- t_mother between (t_lo, t_hi) 
#--------------------------------------------
#--- actual data
#--- remove pairs with r > cs.t
nij, nij_removed = RemovePair( nij, cs, swarm_lohi )
myTlohi = Tlohi(swarm_lohi['date'].min()+datetime.timedelta(days=dt_cut), #--- set cut-off time
               swarm_lohi['date'].max()+datetime.timedelta(days=-dt_cut)
               )
nij = GetTrigListLoHi( nij, 'date',  #--- t_mother > t_lo
                       ( myTlohi.lo, myTlohi.hi ), 
                           swarm_lohi)
nij.reset_index(inplace=True, drop=True)

#--- random data
#for i in xrange(nsamp):
#    nij_rnd = RemovePair( nij_rnd, cs, swarm_shuffled )
#    nij_rnd = GetTrigListLoHi( nij_rnd, 'date', 
#                       ( myTlohi.lo, myTlohi.hi ), 
#                           swarm_shuffled)
#    nij_rnd.reset_index(inplace=True, drop=True)

#--- compute threshold
N_thresh = []
frames = []
for i in xrange(nsamp):
    frames.append(nij_rnd[i])
    #--- plot histograms for random data
    edges, hist, cdf = HIST( nij_rnd[i], 1024 )
#    edges = 0.5*(edges[:-1]+edges[1:])
    N_thresh.append(edges[:-1][cdf>quantile][0])
n_thresh = np.mean(N_thresh)
print 'n_th=%2.1e+%2.1e' %(n_thresh,np.std(N_thresh) / nsamp ** 0.5)
Nij_rnd = pd.concat([nij_rnd[i]]) #frames)


#--- distribution of nij (actual data)
edges_actual, hist_actual, cdf_actual = HIST( nij, 1024 )
#edges_actual = 0.5*(edges_actual[:-1]+edges_actual[1:])
#PLOT(edges, hist, cdf, 'random')

#--- (Tmin,Tmax)-(Rmin,Rmax)
TMIN=np.floor(np.log10(nij['T*']).min())
TMAX=np.ceil(np.log10(nij['T*']).max())
RMIN=np.floor(np.log10(nij['R*']).min())
RMAX=np.ceil(np.log10(nij['R*']).max())

#--------------------
#--- plot scattered
#--------------------
fig = plt.figure(figsize=(8,4))


#--------------------
#--- interpolated data
#--------------------
for nn, index, title in zip([ nij, Nij_rnd ],xrange(2),['actual','randomized']):
#for nn, index, title in zip([ nij ],xrange(2),['actual']):
    ax = fig.add_subplot(1,2,index+1)
    #--- interpolate
    heatmap=GetInterpolatedData(nn,(TMIN,TMAX),(RMIN,RMAX))
    ax.pcolormesh(heatmap.T,cmap='coolwarm')
    #---
#    ax.set_xlabel('T*')
#    ax.set_ylabel('R*')
#    ax.set_title(title)
    ##--- reset x-ticks
    
    reset_xticks(4,ax,(TMIN,TMAX))
    #--- reset y-ticks
    reset_yticks(4,ax,(RMIN,RMAX))
    #--- nij must be within limits
    nij_copy = nn.copy()
    nij_copy=nij_copy[(nij_copy['T*'] >= 10**TMIN ) &
    (nij_copy['T*'] < 10**TMAX ) &
    (nij_copy['R*'] >= 10**RMIN ) &
    (nij_copy['R*'] < 10**RMAX )]
    #--- mapping
    xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
    yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(nij_copy['R*'])-RMIN)/(RMAX-RMIN)
    ax.set_xlim(ax.axis()[0],ax.axis()[1])
    ax.set_ylim(ax.axis()[2],ax.axis()[3])
    #--- plot                 
    DrawFrame(ax, (0.25,0.08),(0.16,0.08),0.01)
    ax.scatter(xx,
              yy,
               alpha=1,color='black',s=0.01)
    #--- draw the separation line
    xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
    yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(n_thresh/nij_copy['T*'])-RMIN)/(RMAX-RMIN)
#    print len(xx)
#    xx=xx[0:-1:300]
#    yy=yy[0:-1:300]
#    print len(xx)
    ax.plot( xx, yy,'-.',color='white',linewidth=1.0)
    
#fig.savefig('%s/densityMaps.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)
#---
# mpl.rcParams.update(mpl.rcParamsDefault)
# warnings.filterwarnings('ignore') #--- get rid of warnings
# rc('text', usetex=True)
# font = {'size'   : 20}
# matplotlib.rc('font', **font)
plt.show()
n_th=2.9e-03+0.0e+00
In [37]:
#--- triggering part of nij
nij_trig=nij[nij['T*']*nij['R*']<=n_thresh]
df_triggered = swarm_lohi.iloc[nij_trig['Event_id']]

#--- background part of nij
nij_background=nij[nij['T*']*nij['R*']>n_thresh]
nn = len(nij_background)
nij_background = nij_background.append(nij_removed,verify_integrity=False) #--- append removed part
assert len(nij_background) == nn + len(nij_removed)

#--- background events
background_child = list(set(nij_background['Event_id'])) #--- children
#
nij_copy = nij.copy()
Nij = nij_copy.append(nij_removed)
background_root  = list(set(Nij['Parent_id'])-set(Nij['Event_id'])) #--- root
#
background = background_root+background_child #--- total
assert len(background) == len(set(background)) #--- assert no overlap
#
df_background = swarm_lohi.iloc[background]

assert len(df_triggered) + len(df_background) == len(swarm_lohi), '%s != %s' %(len(df_triggered) + len(df_background),len(swarm_lohi))


#display(nij_trig.info())
trig_ratio = 1.0*len(df_triggered) / len(swarm_lohi) #cdf_actual[edges_actual[:-1]<=n_thresh][-1]
print 'trig_ratio=%s'%trig_ratio
trig_ratio=0.467287488061

Distributions of $n_{ij}$

In [38]:
#--- histogram
nn=np.log10(nij['R*']*nij['T*'])
nmin=nn.min()
nmax=nn.max()
bins=np.linspace(nmin,nmax,int(nmax-nmin)*8)
hist, edge = np.histogram(nn,bins=bins,normed=True)

#--- accumulated histogram
slist=np.array(nn)
slist.sort()
N = len(slist)
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()

xx=[];yy=[]
for ikey in keys:
    xx.append(d[ikey][0])
    yy.append(d[ikey][2])

#------------
#--- plot
#------------    
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlabel(r'$log n_{ij}$')
ax.set_ylabel(r'$P(log n_{ij})$')
(nmin,nmax) = (min(slist),max(slist))
ax.set_xlim(floor(nmin),ceil(nmax))
#ax.set_xlim(-7,1)
#ax.set_ylim(0,0.6)
#ax.set_xscale('log')
#ax.set_yscale('log')
ax.plot(edge[:-1],hist,'-o',color='black')

# ax2 = ax.twinx()
# ax2.plot(xx,yy,
#         linestyle='-', drawstyle='steps-post',color='red',
#         linewidth=1.0) #--- accumulated
# #ax2.set_xlim(-7,1)
# #ax2.set_ylim(0,1200)

# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')

np.save('actual',np.c_[edge[:-1],hist])
np.save('acc_actual',np.c_[xx,yy])


#fig.savefig('%s/pn_actual.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()

Distributions of random 𝑛𝑖𝑗

In [39]:
#--- histogram
nn=np.log10(Nij_rnd['R*']*Nij_rnd['T*'])
nmin=nn.min()
nmax=nn.max()
bins=np.linspace(nmin,nmax,int(nmax-nmin)*8)
hist, edge = np.histogram(nn,bins=bins,normed=True)

#--- accumulated histogram
slist=np.array(nn)
slist.sort()
N = len(slist)
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()

xx=[];yy=[]
for ikey in keys:
    xx.append(d[ikey][0])
    yy.append(d[ikey][2])

#------------
#--- plot
#------------    
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlabel(r'$log n_{ij}$')
ax.set_ylabel(r'$P(log n_{ij})$')
(nmin,nmax) = (min(slist),max(slist))
ax.set_xlim(floor(nmin),ceil(nmax))
#ax.set_xlim(-7,1)
#ax.set_ylim(0,0.6)
#ax.set_xscale('log')
#ax.set_yscale('log')
ax.plot(edge[:-1],hist,'-o',color='black')

# ax2 = ax.twinx()
# ax2.plot(xx,yy,
#         linestyle='-', drawstyle='steps-post',color='red',
#         linewidth=1.0) #--- accumulated
# #ax2.set_xlim(-7,1)
# #ax2.set_ylim(0,1200)

# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')

np.save('random',np.c_[edge[:-1],hist])
np.save('acc_random',np.c_[xx,yy])


#fig.savefig('%s/pn_shuffled.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()
In [40]:
act = np.load('actual.npy')
rnd = np.load('random.npy')
acc_act = np.load('acc_actual.npy')
acc_rnd = np.load('acc_random.npy')

#------------
#--- plot
#------------    
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlabel(r'$log n_{ij}$')
ax.set_ylabel(r'$P(log n_{ij})$')
(nmin,nmax) = (min(slist),max(slist))
ax.set_xlim(floor(nmin),ceil(nmax))
#ax.set_title('HypoDD')
# ax.set_xlim(-7,1)
# ax.set_ylim(0,0.6)
#ax.set_xscale('log')
#ax.set_yscale('log')
ax.plot(act[:,0],act[:,1],'-o',color='black',label='Actual')
ax.plot(rnd[:,0],rnd[:,1],'-s',color='red',label='Shuffled')
#ax.plot([np.log10(6e-5),np.log10(6e-5)],[0,0.3])
ax.legend(fontsize='xx-small')

# ax2 = ax.twinx()
# ax2.plot(acc_act[:,0],acc_act[:,1],
#         linestyle='-', drawstyle='steps-post',color='black',
#         linewidth=1.0) #--- accumulated
# ax2.plot(acc_rnd[:,0],acc_rnd[:,1],
#         linestyle='-', drawstyle='steps-post',color='red',
#         linewidth=1.0) #--- accumulated
# # ax2.set_xlim(-7,1)
# # ax2.set_ylim(0,1200)

# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')


#fig.savefig('%s/pn_actual_shuffled.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=75)

plt.show()

plot moment tensor and triggered events

In [41]:
# from math import *
# import numpy as np
# import pandas as pd

# def getN(phi_az,alpha):
#     phi_az *= pi/180
#     alpha *= pi/180
    
# #    return np.array([cos(alpha)*cos(phi_az), cos(alpha)*sin(phi_az), sin(alpha) ])
#     return np.array([cos(alpha)*sin(phi_az), cos(alpha)*cos(phi_az), sin(alpha) ]) #--- azimuth with respect to north

# def tensProd(n1,n2):
#     n=3
#     return np.array([[n1[i]*n2[j] for j in xrange(3)] for i in xrange(3)])

# #nij_trig = nij_background #kam

# #trigger_list = set(swarm_lohi.iloc[nij_trig['Parent_id']]['Event ID']) & set(df_moment['Event ID'])
# trigger_list = set(swarm_lohi.iloc[nij_trig['Parent_id']]['Event ID'])


# for triggerID, xx in zip(trigger_list,xrange(10000)):
#     try:
#         #--- list of children
#         index_trigger = swarm_lohi[swarm_lohi['Event ID']==triggerID].index[0]
#         children_index_list = nij_trig[nij_trig['Parent_id']==index_trigger]['Event_id']
#         if len(children_index_list) < 30: 
#             continue
#         #--- plot moment tensor
#         #         #--- parse moment tensor
# #        df_trigger = df_moment[df_moment['Event ID'] == triggerID]
#         df_trigger = swarm_lohi[swarm_lohi['Event ID'] == triggerID]
        
#         n1 = getN(df_trigger['Azimuth(T)'],df_trigger['Plunge(T)']) #--- principal directions (azimuth, plunge)
#         n3 = getN(df_trigger['Azimuth(P)'],df_trigger['Plunge(P)'])
#         n2 = np.cross(n3,n1) #getN(df_trigger['Azimuth(N)'],df_trigger['Plunge(N)'])
#         if np.isnan(n1[0]):
#             continue
#         #--- project t, p
#         n1 -= n1[2] * np.array([0.0,0.0,1.0]) 
#         n3 -= n3[2] * np.array([0.0,0.0,1.0]) 
        
#         #--- plot events
#         fig = plt.figure(figsize=(4,4))
#         ax = fig.add_subplot(111)

#         x = swarm_lohi.iloc[children_index_list]['longitude']
#         y = swarm_lohi.iloc[children_index_list]['latitude']
#         mag = swarm_lohi.iloc[children_index_list]['magnitude']
#         xc = swarm_lohi.iloc[index_trigger]['longitude']
#         yc = swarm_lohi.iloc[index_trigger]['latitude']
#         mag_c = swarm_lohi.iloc[index_trigger]['magnitude']
#         ax.scatter(x,y,s=np.exp(mag),alpha=0.2)
#         ax.scatter(xc,yc,s=np.exp(mag_c),marker='*',facecolors='yellow',color='black')
        
#         dd=max(x.max()-x.min(),y.max()-y.min())
#         rr=.3*dd
#         [x1,y1]=np.array([xc,yc])+rr*n3[:2]
#         [xo,yo]=np.array([xc,yc])+.1*rr*n3[:2]
#         ax.annotate("", (xo,yo), xytext=(x1, y1),
#                          textcoords='data',
#                          arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="red",linewidth="0.3")) 
#         [x1,y1]=np.array([xc,yc])-rr*n3[:2]
#         [xo,yo]=np.array([xc,yc])-.1*rr*n3[:2]
#         ax.annotate("", (xo,yo), xytext=(x1, y1),
#                          textcoords='data',
#                          arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="red",linewidth="0.3")) 
    
#         [x1,y1]=np.array([xc,yc])+rr*n1[:2]
#         [xo,yo]=np.array([xc,yc])+.1*rr*n1[:2]
#         ax.annotate("", (x1,y1), xytext=(xo, yo),
#                          textcoords='data',
#                          arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="black",linewidth="0.3")) 
#         [x1,y1]=np.array([xc,yc])-rr*n1[:2]
#         [xo,yo]=np.array([xc,yc])-.1*rr*n1[:2]
#         ax.annotate("", (x1,y1), xytext=(xo, yo),
#                          textcoords='data',
#                          arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="black",linewidth="0.3")) 
    
        
#         ax.set_xlim(xc-dd/2,xc+dd/2)    
#         ax.set_ylim(yc-dd/2,yc+dd/2)    
        
#         fig.savefig('%s/spatialMap%s.png'%(DIR_OUTPT_figs,triggerID),bbox_inches='tight',dpi=150)
#     except:
#         continue
# #    print (xc,yc),(x1, y1)

# plt.show()
In [42]:
# #trigger_list = set(swarm_lohi.iloc[nij_trig['Parent_id']]['Event ID']) & set(df_moment['Event ID'])


# fig = plt.figure(figsize=(4,4))
# ax = fig.add_subplot(111)
# for triggerID in trigger_list:
#     #--- parse moment tensor
# #    df_trigger = df_moment[df_moment['Event ID'] == triggerID]
#     df_trigger = swarm_lohi[swarm_lohi['Event ID'] == triggerID]
        
#     n1 = getN(df_trigger['Azimuth(T)'],df_trigger['Plunge(T)']) #--- principal directions (azimuth, plunge)
#     n3 = getN(df_trigger['Azimuth(P)'],df_trigger['Plunge(P)'])
#     n2 = np.cross(n3,n1) #getN(df_trigger['Azimuth(N)'],df_trigger['Plunge(N)'])
#     if np.isnan(n1[0]):
#             continue
            
#     Q=np.c_[n1,n3,n2].T #--- transformation 
    
#     #--- list of children
#     index_trigger = swarm_lohi[swarm_lohi['Event ID']==triggerID].index[0]
#     children_index_list = nij_trig[nij_trig['Parent_id']==index_trigger]['Event_id']
#     #--- plot events
    
#     x = swarm_lohi.iloc[children_index_list]['x(km)']
#     y = swarm_lohi.iloc[children_index_list]['y(km)']
#     z = swarm_lohi.iloc[children_index_list]['z(km)']
#     xc = swarm_lohi.iloc[index_trigger]['x(km)']
#     yc = swarm_lohi.iloc[index_trigger]['y(km)']
#     zc = swarm_lohi.iloc[index_trigger]['z(km)']
    
#     r_local = np.matmul(Q,np.c_[x-xc,y-yc,z-zc].T).T #Transform
    
#     ax.scatter(r_local[:,0],r_local[:,1],c='C0',alpha=0.1)
# #    ax.scatter(xc,yc,color='yellow')
#     #--- plot moment tensor

    
#     #--- project t, p
#     n1 = np.array([1.0,0.0]) 
#     n3 = np.array([0.0,1.0]) 
#     rr=4 #.2*(r_local[:,0].max()-r_local[:,0].min())
    
# xc = 0.0
# yc = 0.0
# zc = 0.0
# [x1,y1]=np.array([xc,yc])+rr*n3[:2]
# ax.annotate("", (xc,yc), xytext=(x1, y1),
#                  textcoords='data',
#                  arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="red",linewidth="0.3")) 
# [x1,y1]=np.array([xc,yc])-rr*n3[:2]
# ax.annotate("", (xc,yc), xytext=(x1, y1),
#                  textcoords='data',
#                  arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="red",linewidth="0.3")) 

# [x1,y1]=np.array([xc,yc])+rr*n1[:2]
# ax.annotate("", (x1,y1), xytext=(xc, yc),
#                  textcoords='data',
#                  arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="black",linewidth="0.3")) 
# [x1,y1]=np.array([xc,yc])-rr*n1[:2]
# ax.annotate("", (x1,y1), xytext=(xc, yc),
#                  textcoords='data',
#                  arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.4",color="black",linewidth="0.3")) 

# dd=20 #max(r_local[:,0].max()-r_local[:,0].min(),r_local[:,1].max()-r_local[:,1].min())
# ax.set_xlim(xc-dd/2,xc+dd/2)    
# ax.set_ylim(yc-dd/2,yc+dd/2)    
# ax.set_xlabel(r'$x(km)$')
# ax.set_ylabel(r'$y(km)$')
# fig.savefig('%s/spatialMap.png'%(DIR_OUTPT_figs),bbox_inches='tight',dpi=150)

# plt.show()

Declustered Catalog

In [43]:
df_triggered.head()
Out[43]:
date latitude longitude depth magnitude ID x(km) y(km) z(km)
1 2019-07-04 17:09:20.000 35.709648 -117.499463 12.435 2.49 70461760 35.414181 26.034539 12.435
2 2019-07-04 17:12:14.700 35.708114 -117.499284 11.876 2.13 70461935 35.434107 25.864336 11.876
4 2019-07-04 17:35:01.640 35.643070 -117.557894 10.229 4.49 70463302 28.909672 18.647509 10.229
5 2019-07-04 17:35:52.000 35.687642 -117.509147 0.830 3.00 70463352 34.336163 23.592898 0.830
6 2019-07-04 17:36:38.520 35.723421 -117.537760 3.545 2.80 70463397 31.150978 27.562705 3.545
In [44]:
#--- plot complete catalog

def ConvertDailyRate(hist, bin_edges ):
#---convert data to daily rate     
    t0 = pd.to_datetime( bin_edges[ 0 ] )
    t1 = pd.to_datetime( bin_edges[ 1 ] )
    delta_t = ( t1 - t0 ).total_seconds() / ( 60 * 60 * 24)
    hist *= ( bin_edges[ 1 ] - bin_edges[ 0 ] ) / delta_t

def ActivityRate( df ):
    nbins = int( (df['date'].max()-df['date'].min()).days / 2 ) #--- number of bins
    
    tmax = df['date'].max().value #--- min/max
    tmin = df['date'].min().value
    
    hist, bin_edges = np.histogram(df['date'].apply(lambda x:x.value),                                   
                    bins=np.linspace(tmin,tmax,nbins+1,endpoint=True),density=True) #--- probability dist.
    hist *= len( df['date'] ) #--- int(hist).dt=n
    cumm_number = np.cumsum(hist)*(bin_edges[1]-bin_edges[0]) #--- accumulated number
    ConvertDailyRate( hist, bin_edges ) #--- daily activity
    return bin_edges, hist, cumm_number
    
class Tlohi:
    def __init__(self,lo,hi):
        self.lo = lo
        self.hi = hi
        
#---------------------------------------------------------------------------------
#-----------------
#-----------------
#-----------------
#---------------------------------------------------------------------------------
#--- set min/max time to avoid temporal incompletenesss issue
df_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
df_lohi = GetCartesian( df_lohi )
df_lohi.reset_index(inplace=True,drop=True)

#--- activity rate
bin_edges, hist, cumm_number = ActivityRate( df_lohi )

#--- temporal map
fig = plt.figure(figsize=(8,9))#,dpi=150)
ax1=fig.add_subplot(311)
ax2=fig.add_subplot(312)
ax3=fig.add_subplot(313)
ax1.set_ylabel('$M$')
ax1.set_xlim(df_lohi['date'].min(),df_lohi['date'].max())
ax1.tick_params(axis='both',labelsize=18)
ax2.set_ylabel('$M$')
ax2.set_xlim(df_lohi['date'].min(),df_lohi['date'].max())
ax2.tick_params(axis='both',labelsize=18)
ax3.set_ylabel('$M$')
ax3.set_xlim(df_lohi['date'].min(),df_lohi['date'].max())
ax3.tick_params(axis='both',labelsize=18)
#ax1.set_ylim(mc,2.5)
#ax1.set_xlim(pd.to_datetime('2019-03-01'),pd.to_datetime('2019-03-30'))
#
# ax1R = ax1.twinx()  # instantiate a second axes that shares the same x-axis
# ax1R.set_ylabel('daily activity rate',color='red')
# ax1R.set_yscale('log')
# ax1R.tick_params(axis='y',labelsize=18,colors='red')
#
#ax2 = fig.add_subplot(212)
#ax2.set_ylim(mc,2.5)
#ax2.set_ylabel('Max Rate (m3/min)')
#ax2.set_ylim(11.9,12.1)
#ax2.set_xlim(pd.to_datetime('2019-03-01'),pd.to_datetime('2019-03-30'))
#

#--- plot
ax1.scatter(df_triggered['date'],df_triggered['magnitude'],
            s=2**(df_triggered['magnitude']),
            alpha=0.05,color='blue', label = 'Triggered')
#
ax2.scatter(df_background['date'],df_background['magnitude'],
            s=2**(df_background['magnitude']),
            alpha=0.05,color='red', label = 'Declustered')
#
ax3.scatter(df_lohi['date'],df_lohi['magnitude'],
            s=2**(df_lohi['magnitude']),
            alpha=0.05,color='black', label = 'Total')
fig.legend(loc='best',frameon=True,labelspacing=0)
#ax1R.plot(pd.to_datetime(bin_edges[:-1]),hist,'r-',linewidth=1)

#DrawFrame(ax1, (0.11,0.07),(0.32,0.1),0.01)
ax3.xaxis.set_major_formatter(dates.DateFormatter('%b %y'))
fig.autofmt_xdate()

fig.savefig('%s/timeSeries_ok.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')

plt.show()

TypeErrorTraceback (most recent call last)
<ipython-input-44-de18f5f41f06> in <module>()
     73 ax1.scatter(df_triggered['date'],df_triggered['magnitude'],
     74             s=2**(df_triggered['magnitude']),
---> 75             alpha=0.05,color='blue', label = 'Triggered')
     76 #
     77 ax2.scatter(df_background['date'],df_background['magnitude'],

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1715                     warnings.warn(msg % (label_namer, func.__name__),
   1716                                   RuntimeWarning, stacklevel=2)
-> 1717             return func(ax, *args, **kwargs)
   1718         pre_doc = inner.__doc__
   1719         if pre_doc is None:

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, **kwargs)
   4021             linewidths = rcParams['lines.linewidth']
   4022 
-> 4023         offsets = np.column_stack([x, y])
   4024 
   4025         collection = mcoll.PathCollection(

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/numpy/lib/shape_base.pyc in column_stack(tup)
    367             arr = array(arr, copy=False, subok=True, ndmin=2).T
    368         arrays.append(arr)
--> 369     return _nx.concatenate(arrays, 1)
    370 
    371 def dstack(tup):

TypeError: invalid type promotion
In [45]:
#--- plot spatial map 
#--- fault map california: temporal evolution of events


#--- set min/max time to avoid temporal incompletenesss issue
df_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
df_lohi = GetCartesian( df_lohi )
df_lohi.reset_index(inplace=True,drop=True)
# df_triggered = GetCartesian( df_triggered )
# df_background = GetCartesian( df_background )


xlo = df_lohi['x(km)'].min()
xhi = df_lohi['x(km)'].max()
ylo = df_lohi['y(km)'].min()
yhi = df_lohi['y(km)'].max()
L = max(xhi-xlo,yhi-ylo)
exp = 2

df_triggered.plot.scatter('x(km)','y(km)',
                        s=exp**(df_triggered['magnitude']),
                        color = 'blue', label = 'triggered',xlim=(0,L),ylim=(0,L),
#                        c='date',cmap='jet',
                        alpha=0.01,
                        figsize=(4,4)) #--- plot

df_background.plot.scatter('x(km)','y(km)',
                        s=exp**(df_background['magnitude']),
                        color = 'red', label = 'background',xlim=(0,L),ylim=(0,L),
#                        c='date',cmap='jet',
                        alpha=0.01,
                        figsize=(4,4)) #--- plot

df_lohi.plot.scatter('x(km)','y(km)',
                        s=exp**(df_lohi['magnitude']),
                        color = 'black', label = 'Total',xlim=(0,L),ylim=(0,L),
#                        c='date',cmap='jet',
                        alpha=0.01,
                        figsize=(4,4)) #--- plot







plt.show()
plt.savefig('%s/map.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
In [46]:
#trig_ratio *= 0.75
#n_thresh_new = edges_actual[:-1][cdf_actual <= trig_ratio][-1]
#new_quantile = cdf[edges[:-1]<=n_thresh_new][-1]
#print 'new_quantile=%s'%new_quantile
#n_thresh = n_thresh_new
In [47]:
#--- plot clusters


def Inside(t,(tmin,tmax)):
    if tmin<= t<tmax:
        return True
    return False    

def GetTrigList( nij_trig ):
    d_trig = {}
    for items in nij_trig.itertuples():
        triggerID = items.Parent_id
        d_trig.setdefault( triggerID, [] ).append( items.Event_id )
    return d_trig

def PlotArrowsTime(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
    tlist=[]
    mlist=[]
    for triggerID in d_trig: #--- mother events
        
        t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
        m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
        if not ( Inside(t0,(tmin,tmax) ) and #--- within the given interval
                 Inside(swarm_lohi['x(km)'].iloc[triggerID],(xmin,xmax)) and 
                 Inside(swarm_lohi['y(km)'].iloc[triggerID],(ymin,ymax)) ):
            continue

#        x0 = md.date2num( t0 ) #--- convert to time object
#        y0 = m0 #--- magnitude
        
        tlist.append( t0 )
        mlist.append( m0 )

        for daughter_id in d_trig[triggerID]: #--- children
            
            t1 = df_complete['date'].iloc[ daughter_id ] #--- time
            m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
            
            tlist.append( t1 )
            mlist.append( m1 )
            
            if not ( Inside(t1,(tmin,tmax) ) and #--- within the given interval
                     Inside(swarm_lohi['x(km)'].iloc[daughter_id],(xmin,xmax)) and 
                     Inside(swarm_lohi['y(km)'].iloc[daughter_id],(ymin,ymax)) ):
                continue
                
#            x1 = md.date2num( t1 ) #--- convert to time object
#            y1 = m1 #--- magnitude

#            xw = x1 - x0 #--- t_child-t_parent
#            yw = y1 - y0 #--- m_child - m_parent
            plt.annotate("", (t1,m1), xytext=(t0, m0),
                         textcoords='data',
                        arrowprops=dict(arrowstyle="-|>,head_width=.4,head_length=.8",color="b",linewidth="0.3")) 
    
    #--- plot circles
    df=pd.DataFrame({'date':tlist,'mag':mlist})
    plt.scatter(df['date'], df['mag'],
                s=4**(df['mag']),
                alpha=1,
                facecolors='red',color='black')
    #plt.savefig('timeSeries.png')
    
#--------------------------------------
#----- key: event value: aftershock id
#--------------------------------------
d_trig = GetTrigList( nij_trig )

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)

#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )



# #--- plot arrows
# #dt=30
# for ii in xrange(1):
#     #--- setup
#     fig = plt.figure(figsize=(16,4),frameon=False) #, dpi=300)
#     #--- xlimit
#     xmin = swarm_lohi['longitude'].min()
#     xmax = swarm_lohi['longitude'].max()
#     ymin = swarm_lohi['latitude'].min()
#     ymax = swarm_lohi['latitude'].max()
#     plt.tick_params(axis='x', labelsize=14)
#     plt.tick_params(axis='y', labelsize=14)
    
#     #--- plot all events
#     plt.ylim(mc,swarm_lohi['magnitude'].max())
#     plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
#                 s=2*np.exp(swarm_lohi['magnitude']),
#             alpha=0.1)

#     key_max = swarm_lohi.loc[nij_trig.groupby(by='Parent_id').groups.keys()]['magnitude'].sort_values(ascending=False).index[ii]
# #    print 'key=',key_max
# #    print 't=',swarm_lohi.loc[key_max]['date']
# #    print 'x=',swarm_lohi.loc[key_max]['x(km)'],
# #    print 'y=',swarm_lohi.loc[key_max]['y(km)'],
# #    print 'z=',swarm_lohi.loc[key_max]['z(km)']
#     dt = nij_trig[nij_trig['Parent_id']==key_max]['T'].max()
#     dx = nij_trig[nij_trig['Parent_id']==key_max]['R'].max()
#     tt=swarm_lohi['date'].loc[key_max]
#     xx = swarm_lohi.loc[key_max][['x(km)','y(km)','z(km)']]
#     PlotArrowsTime( swarm_lohi, d_trig,
#                        ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)), 
#            (xx['x(km)']-dx,xx['x(km)']+dx),
#            (xx['y(km)']-dx,xx['y(km)']+dx) )
#     plt.xlim(tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt))
In [48]:
# #--- plot clusters


# def Inside(t,(tmin,tmax)):
#     if tmin<= t<tmax:
#         return True
#     return False    

# def GetTrigList( nij_trig ):
#     d_trig = {}
#     for items in nij_trig.itertuples():
#         triggerID = items.Parent_id
#         d_trig.setdefault( triggerID, [] ).append( items.Event_id )
#     return d_trig

    
# def PlotArrowsSpace(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
#     xlist=[]
#     ylist=[]
#     mlist=[]
#     for triggerID in d_trig: #--- mother events
        
#         t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
#         ( x0, y0 ) = ( swarm_lohi['longitude'].iloc[triggerID], 
#                        swarm_lohi['latitude'].iloc[triggerID] )
#         m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
# #        if not ( Inside(t0,(tmin,tmax) ) and #--- within the given interval
#  #                Inside(x0,(xmin,xmax)) and 
#   #               Inside(y0,(ymin,ymax)) ):
#    #         continue

#         for daughter_id in d_trig[triggerID]: #--- children
#             t1 = df_complete['date'].iloc[ daughter_id ] #--- time
#             ( x1, y1 ) = ( swarm_lohi['longitude'].iloc[ daughter_id ], 
#                            swarm_lohi['latitude'].iloc[ daughter_id ] )
#             m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
            
#             xlist.append( x1 )
#             ylist.append( y1 )
#             mlist.append( m1 )
            
# #            if not ( Inside(t1,(tmin,tmax) ) and #--- within the given interval
# #                     Inside(x1,(xmin,xmax)) and 
# #                     Inside(y1,(ymin,ymax)) ):
# #                continue

#             ax.annotate("", (x1,y1), xytext=(x0, y0),
#                          textcoords='data',
#                         arrowprops=dict(arrowstyle="-|>,head_width=.2,head_length=.8",color="b",linewidth="0.3")) 
 
#         #--- plot mother
#         plt.scatter(x0,y0,
#                     s=exponent**m0,
#                     marker='*',
#                     facecolors='yellow',color='black',
#                     alpha=1.0)

#     #--- plot circles
#     df=pd.DataFrame({'x(km)':xlist,'y(km)':ylist,'m':mlist})
#     ax.scatter(df['x(km)'], df['y(km)'],
#                 s=exponent**(df['m']),
#                 alpha=1,
#                 facecolors='red',color='black')
    
# #--------------------------------------
# #----- key: event value: aftershock id
# #--------------------------------------
# d_trig = GetTrigList( nij_trig )

# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# swarm_lohi.reset_index(inplace=True,drop=True)

# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )



# #--- plot arrows
# exponent = 4
# for ii in xrange(1):
#     #--- setup
#     fig = plt.figure(figsize=(8,8)) #,frameon=False) #, dpi=300)
#     ax = fig.add_subplot(111)
#     #--- xlimit
#     xmin = swarm_lohi['longitude'].min()
#     xmax = swarm_lohi['longitude'].max()
#     ymin = swarm_lohi['latitude'].min()
#     ymax = swarm_lohi['latitude'].max()
#     ax.tick_params(axis='x', labelsize=14)
#     ax.tick_params(axis='y', labelsize=14)
    
#     #--- plot all events
#     ax.scatter(swarm_lohi['longitude'], swarm_lohi['latitude'], # 'x(km)'],swarm_lohi['y(km)'],
#                 s=exponent**swarm_lohi['magnitude'],
#                 alpha=0.1)

# #    key_max = swarm_lohi.loc[nij_trig.groupby(by='Parent_id').groups.keys()]['magnitude'].sort_values(ascending=False).index[ii]
#     key_max = nij_trig.groupby(by='Parent_id').count()['R'].sort_values(ascending=False).index[ii]
#     dt = nij_trig[nij_trig['Parent_id']==key_max]['T'].max()
#     dx = 0.1 #nij_trig[nij_trig['Parent_id']==key_max]['R'].max()
#     tt = swarm_lohi['date'].loc[key_max]
#     xx = swarm_lohi.loc[ key_max ][['longitude','latitude']] #[['x(km)','y(km)','z(km)']]
    
#     PlotArrowsSpace( swarm_lohi, d_trig,
#                        ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)), 
#            (xx['longitude']-dx,xx['longitude']+dx),
#            (xx['latitude']-dx,xx['latitude']+dx) )
#     ax.set_xlim(xx[0]-dx,xx[0]+dx)
#     ax.set_ylim(xx[1]-dx,xx[1]+dx)
In [49]:
# #--- spatial clusters



# #def GetCartesian2nd( dff ):
# #    df = dff.copy()
# #    xlo = df['longitude'].min()
# #    xhi = df['longitude'].max()
# #    ylo = df['latitude'].min()
# #    yhi = df['latitude'].max()
# #    getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
# #    getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
# #    df[ 'x' ] = df[ 'longitude' ].apply( getDistX )
# #    df[ 'y' ] = df[ 'latitude' ].apply( getDistY )
# #    return df


# def PlotMap(df_complete,TRIG_ID, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
#     #--- plot 
#     counter = 0
#     exponent=4
#     for triggerID,colin in zip(TRIG_ID,colors):
#         t0 = df_complete['date'].iloc[ triggerID ]
# #        if not ( Inside(t0,(tmin,tmax) ) and 
# #                 Inside(df_complete['y(m)'].iloc[triggerID],(xmin,xmax)) and 
# #                 Inside(df_complete['z(m)'].iloc[triggerID],(ymin,ymax)) ):  
# #              continue
#         counter += 1
#         if counter > 10:
#             break
#         #--- plot    
#         fig = plt.figure(figsize=(5,5))#,dpi=150)
#         ax = fig.add_subplot(111)
# #        plt.xlim(xmin,xmax)
# #        plt.ylim(ymin,ymax)
#         plt.tick_params(axis='x', labelsize=14)
#         plt.tick_params(axis='y', labelsize=14)
#         plt.scatter(df_complete.loc[triggerID]['x(km)'],
#                     df_complete.loc[triggerID]['y(km)'],
#                     s=exponent**df_complete.loc[triggerID]['magnitude'],
#                     facecolors='white',color='white')
#         d={'x':[],'y':[],'m':[],'t':[]}
#         #--- plot daughters
#         for daughter_id in d_trig[triggerID]:
#             t1 = df_complete['date'].iloc[ daughter_id ]
# #            if not ( Inside(t1,(tmin,tmax) ) and 
# #                     Inside(df_complete['y(m)'].iloc[daughter_id],(xmin,xmax)) and 
# #                     Inside(df_complete['z(m)'].iloc[daughter_id],(ymin,ymax)) ):
# #                continue
#     #--- remove pairs with r > cs.t
# #            if rmat[ triggerID, daughter_id ] > cs * tmat[ triggerID, daughter_id ]:
# #                print 'hello!'
# #                continue

#             d['t'].append( t1 )
#             d['x'].append(df_complete['x(km)'].iloc[daughter_id])
#             d['y'].append(df_complete['y(km)'].iloc[daughter_id])
#             d['m'].append(df_complete['magnitude'].iloc[daughter_id])
#         df=pd.DataFrame(d)
#         plt.scatter(df['x'],df['y'],
#                     s=exponent**df['m'],
#                     c=df['t'], cmap='jet') #,alpha=0.1)
#         plt.colorbar()
#         #--- plot mother
#         plt.scatter(df_complete['x(km)'].iloc[triggerID],df_complete['y(km)'].iloc[triggerID],
#                     s=exponent**df_complete['magnitude'].iloc[triggerID],
#                     marker='*',
#                     facecolors='yellow',color='black',
#                     alpha=1.0)
#         #--- 


    
# colors = itertools.cycle(["r", "b", "g","b"])
# #cs = 3 * (24*3600) #--- km per day

# #--- sort
# new_list = [[len(d_trig[triggerID]),triggerID] for triggerID in d_trig]
# new_list.sort(reverse=True)
# TRIG_ID = [i[1] for i in new_list]

# #--- xlimit
# swarm_lohi=GetCartesian(swarm_lohi)
# #swarm_lohi['y(m)']=swarm_lohi['x']
# #swarm_lohi['z(m)']=swarm_lohi['y']
# xmin = swarm_lohi['x(km)'].min()
# xmax = swarm_lohi['x(km)'].max()
# ymin = swarm_lohi['y(km)'].min()
# ymax = swarm_lohi['y(km)'].max()
# #plt.tick_params(axis='x', labelsize=14)
# #plt.tick_params(axis='y', labelsize=14)

# #--- plot arrows
# PlotMap( swarm_lohi, TRIG_ID,
# #          (pd.to_datetime('2010-04-04'),pd.to_datetime('2010-06-26')),
# #          (pd.to_datetime('2010-04-08'),pd.to_datetime('2010-04-14')),
#                        ( swarm['date'].min(), swarm['date'].max()), 
#            (xmin,xmax),
#            (ymin,ymax) )
In [50]:
#--- cluster life time


def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
    df_parents = df.groupby(by='Parent_id').groups #--- parent events
    ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
    ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
    df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
    return df_m

#def AddDist( df_trig, df_complete ):
#    x=df_complete.loc[ df_trig['Event_id'] ]['r(m)'] 
#    y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
#    df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
#    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )

def AddTime( df_trig, df_complete ):
    prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
    x=df_complete.loc[ df_trig['Event_id'] ]
    y=df_complete.loc[ df_trig['Parent_id'] ]
    x.reset_index(inplace=True)
    y.reset_index(inplace=True)
    df_trig['T'] = x['date']-y['date']
    df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)

def RemovePair( df, cs ):
    return df [ df[ 'R' ] <= df[ 'T' ] * cs ]

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)


#--- choose m0<m<m1
n = 7
#m = list(np.linspace(mc,swarm['magnitude'].max(),n))
m_parents = swarm_lohi.loc[nij_trig.groupby(by='Parent_id').groups.keys()]['magnitude']
m = list(np.linspace(m_parents.min()-1e-6,m_parents.max()+1e-6,n))
#m+=[swarm['magnitude'].max()]

X={}
rho_dict={}
m_range={}
Err={}
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
for i in xrange( len(m) - 1 ):
    m0 = m[ i ]
    m1 = m[ i + 1 ]
        
#    print m0, m1
#    nij_trig.groupby(by='Parent_id')
    
    df_trig = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
    df_trig.reset_index( inplace = True, drop = True )
    #--- add distance & time column
#    df_trig['R'] = rmat[df_trig['Parent_id'],df_trig['Event_id']]
#    df_trig['T'] = tmat[df_trig['Parent_id'],df_trig['Event_id']]
#    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
#    assert len ( df_trig[ df_trig['T'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['T'] == 0.0 ] )
    
    #--- remove pairs with r > cs.t
#    cs = 3 * (24*3600) #--- km per day
#    df_trig = RemovePair( df_trig, cs )

    #--- life time
    df_duration = df_trig.groupby(by='Parent_id')['T'].max()

    #--- rho plots
    n_decades = int( math.ceil( log(df_duration.max()/df_duration.min(),10) ) )
    nbin_per_decade = 3
    nbins =  nbin_per_decade * n_decades
    rho, xedges=np.histogram(df_duration, 
                             density=True,
                             bins=np.logspace(log(df_duration.min(),10),
                                              log(df_duration.max(),10),nbins))
    hist, xedges=np.histogram(df_duration,
                              bins=np.logspace(log(df_duration.min(),10),
                                               log(df_duration.max(),10),nbins))
    err = rho/np.sqrt(hist)
        
    
    #--- rho
#    plt.xlim(1e-3,1e2)
#    plt.ylim(1e-3,1e2)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('Duration(Day)')
    plt.ylabel('P(Duration)')
    plt.errorbar(xedges[:-1],rho, yerr=err,label='[%2.1f,%2.1f]'%(m0,m1))
    fontP = FontProperties()
    fontP.set_size('small')
    plt.legend(title='m',prop=fontP,loc='upper right',bbox_to_anchor=(1.5, 1))
    x=xedges[:-1]
    
    X[i] = x[:]
    rho_dict[i]=rho[:]
    m_range[i]=np.array([m0,m1])
    Err[i]=err[:]
#plt.savefig('%s/duration.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
In [51]:
#--- histograms: magnitude of mother events

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)


mlist = swarm_lohi.loc[ nij_trig.groupby(by='Parent_id').groups.keys() ]['magnitude']
#plt.hist(mlist)

list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
mlist0= swarm_lohi['magnitude'].loc[df_zeros.index]


#--- accumulated histogram
N = len(mlist)+len(mlist0)
slist=mlist.tolist()+mlist0.tolist() #np.array(swarm_copy['magnitude'])
slist.sort()
d = histogramACCUMLTD( slist )
keys=d.keys()
keys.sort()

#--- plot
fig, ax = plt.subplots(figsize=(4,4))
ax.set_yscale('log')
ax.set_xlabel('$m$')
ax.set_ylabel('$N(m)$')
xmin=min(slist)
xmax=max(slist)
ax.axis([xmin,xmax,1,1000])
#ax.set_xticks(np.linspace(1,6,6))

dx=0.1
junk = ax.hist( slist,
                bins=int((xmax-xmin)/dx),
               label='histogram',color='silver') #--- histogram
xx=[];yy=[]
for ikey in keys:
    xx.append(d[ikey][0])
    yy.append(N-d[ikey][2])
ax.plot(xx,yy,
            linestyle='-', drawstyle='steps-post',color='black',
             linewidth=1.0) #--- accumulated
    
b=1.65
ax.plot([xmin,xmax],
        [1e3, 1e3* 10**(-b*(xmax-xmin))],'r-.',linewidth=1)

#DrawFrame(ax, (0.2,0.06),(0.14,0.06),0.1,LOG_Y=True)

#plt.savefig('%s/gr.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
Out[51]:
[<matplotlib.lines.Line2D at 0x7fe74a6f3f90>]
In [52]:
#--- productivity relation

mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings

rc('text', usetex=True)
font = {'size'   : 26}
matplotlib.rc('font', **font)

import numpy as np
import math

def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
    df_parents = df.groupby(by='Parent_id').groups #--- parent events
    ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
    ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
    df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
    return df_m, np.mean( ds_m ), np.std( ds_m )/(len(ds_m))**0.5 #--- exponential: std=mean

#def AddDist( df_trig, df_complete ):
#    x=df_complete.loc[ df_trig['Event_id'] ]['r(m)'] 
#    y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
#    df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
#    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )

def AddTime( df_trig, df_complete ):
    prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
    x=df_complete.loc[ df_trig['Event_id'] ]
    y=df_complete.loc[ df_trig['Parent_id'] ]
    x.reset_index(inplace=True)
    y.reset_index(inplace=True)
    df_trig['T'] = x['date']-y['date']
    df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)

def RemovePair( dff, cs ):
    df = dff.copy()
    df['R'] = rmat[df['Parent_id'],df['Event_id']]
    df['T'] = tmat[df['Parent_id'],df['Event_id']]
    assert len ( df[ df['R'] == 0.0 ] ) == 0, '%s'%display( df[ df['R'] == 0.0 ] )
    assert len ( df[ df['T'] == 0.0 ] ) == 0, '%s'%display( df[ df['T'] == 0.0 ] )
    return df [ df[ 'R' ] <= df[ 'T' ] * cs ]

def AdjustBins(data,bins):
    hist, xedges=np.histogram(data, #--- histogram
                              bins=bins)
    tmp = [xedges[0]]
    for i,j in zip(xedges[1:],hist): #--- expand bins 
        if j < 2:
            continue
        tmp.append(i)
#    tmp.append(xedges[-1])
    if hist[-1] < 2:
        tmp.pop(-2)
    return tmp 

def GetParentsMagnitudes( swarm_lohi, df_trig ):
    mlist = swarm_lohi.loc[ df_trig.groupby(by='Parent_id').groups.keys() ]['magnitude']

    list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
    df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
    mlist0= swarm_lohi['magnitude'].loc[df_zeros.index]

    return mlist.tolist()+mlist0.tolist()

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)



#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()

#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )

#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']


#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)


#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
#ax.set_xlim(mc,swarm_lohi['magnitude'].max()+0.5)
m_min = min(m_list)-1.0e-6
m_max = max(m_list)+1.0e-6
ax.set_xlim(np.floor(m_min),np.ceil(m_max))
ax.set_ylim(0.2,10**np.ceil(np.log10(df_nn.max())))
ax.set_yscale('log')
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
           df_nn,
          'x',color='black',alpha=0.25,zorder=1,label='')
ax.legend(fontsize=0)

ax2=ax.twinx()
ax2.set_yscale('linear')
ax2.set_xlim(np.floor(m_min),np.ceil(m_max))
ax2.set_ylim(-.02,1)
ax2.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
         df_zeros,
          'x',color='black',alpha=0.25)
ax2.tick_params(labelright=False,labelleft=False,left=False, right=False)

#plt.yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')

#--- choose m0<m<m1
dm = 0.1
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
m=AdjustBins(m_list,m)
Mlist=m

m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
    try:
        m0 = m[ i ]
        m1 = m[ i + 1 ]
        df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
                                                         swarm_lohi,df_trig) #--- get parent with m0<m<m1
        df.reset_index( inplace = True, drop = True )
        
        #--- number of daughters for each mother
        df_n = df.groupby(by='Parent_id').count()['Event_id']
    
        #--- include daughters with zero children
        list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
        subset = swarm_lohi.loc[list_of_zero_afs]
        list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index    
        df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
        
        N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
        err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
#        err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
    except:
        print traceback.print_exc()
        continue
    
#ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
            yerr=err_N,xerr=m_err,
            markerfacecolor='white',
            fmt='o', markersize=8,color='red',markeredgewidth=0.7,
            linewidth=.5,
            barsabove=None,capsize=5,capthick=1,elinewidth=1,
           zorder=2)


#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))

DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.01,LOG_Y=True)

#from scipy.optimize import curve_fit
def func(x,alpha,c):
    return c*10**(alpha*x)
#popt, pconv = curve_fit(func,m_mean,N_as)

popt = np.polyfit(m_mean,np.log10(N_as), 1)
popt[1] = 10**(popt[1])
alpha_p = popt[0]

mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
     swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
mm=list(set(mm))
mm.sort()
mm=np.array(mm)

ax.plot(mm,func(mm,*popt),'-.g',zorder=3,label='$N\propto 10^{%3.2f m}$'%popt[0])
ax.legend(fontsize='xx-small',frameon=False)
ax.set_title('$q=%d,d_f=%s$'%(quantile*100,Df),fontsize='xx-small')
#print 'alpha=%s'%popt[0]
#print pconv**0.5

plt.savefig('%s/prod.png'%DIR_OUTPT_figs,dpi=75*.75,bbox_inches='tight')

mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings

rc('text', usetex=True)
font = {'size'   : 20 }
matplotlib.rc('font', **font)

np.save('%s/alpha'%DIR_OUTPT_figs,alpha_p)


plt.show()
No handlers could be found for logger "matplotlib.legend"
<matplotlib.figure.Figure at 0x7fe774f04e10>
In [53]:
raise SystemExit("Stop right there!")
An exception has occurred, use %tb to see the full traceback.

SystemExit: Stop right there!

multiple realizations

In [54]:
from IPython.display import Image

alpha_arr=[]
for irun in xrange(16):
    path = 'runFillmore/Run%s'%irun
    a=np.load('%s/alpha.npy'%path)
    alpha_arr.append(a.item())
    #--- load image
    display(Image(filename='%s/prod.png'%path))
df=pd.DataFrame({'mean':[np.mean(alpha_arr)],'std_dev':np.std(alpha_arr),
                        'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5]},
                          index=[r'$\alpha$'])
display(df)
df['mean']+2*df['std_dev'],df['mean']-2*df['std_dev']                       
                       

IOErrorTraceback (most recent call last)
<ipython-input-54-06e92f3df818> in <module>()
      4 for irun in xrange(16):
      5     path = 'runFillmore/Run%s'%irun
----> 6     a=np.load('%s/alpha.npy'%path)
      7     alpha_arr.append(a.item())
      8     #--- load image

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/numpy/lib/npyio.pyc in load(file, mmap_mode, allow_pickle, fix_imports, encoding)
    370     own_fid = False
    371     if isinstance(file, basestring):
--> 372         fid = open(file, "rb")
    373         own_fid = True
    374     elif is_pathlib_path(file):

IOError: [Errno 2] No such file or directory: 'runFillmore/Run0/alpha.npy'
In [55]:
#--- number distributions


#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)



#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()

#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )

#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']


#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)

#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
           df_nn,
          'x',color='gray',alpha=1.0)
ax.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
         df_zeros,
          'x',color='gray')
ax.set_yscale('log')


#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )

#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']


#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)


#--- scattered plot
m_min = min(m_list)-1e-6 #0.05 #-1.0e-6
m_max = max(m_list)+1.0e-6 #3.299999333333333 #max(m_list)+0.05 #+1.0e-6 #--- exclude the largest event

#--- choose m0<m<m1
dm = 0.4
n = int( ( m_max - m_min ) / dm )
m = [m_min,m_max] #list(np.linspace(m_min,m_max,n+1))
Mlist=m

m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange(len(m) - 1 ):
    try:
        m0 = m[ i ]
        m1 = m[ i + 1 ]
        df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
                                                         swarm_lohi,df_trig) #--- get parent with m0<m<m1
        df.reset_index( inplace = True, drop = True )
        
        #--- number of daughters for each mother
        df_n = df.groupby(by='Parent_id').count()['Event_id']
    
        #--- include daughters with zero children
        list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
        subset = swarm_lohi.loc[list_of_zero_afs]
        list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index    
        df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
        
#        fig.hist(df_n)
#        if m1 > 3.3:
#            continue
        fig = plt.figure(figsize=(4,4),frameon=True)
        ax = fig.add_subplot(111)
        junk = ax.hist( df_n.tolist()+df_zeros.tolist(),
                bins=1024,#int((xmax-xmin)/dx),
               label='%3.2f-%3.2f'%(m0,m1),color='silver') #--- histogram
        print '#=',junk[0][junk[0]!=0]
        print 'N=',junk[1][:-1][junk[0]!=0]
        ax.set_xlabel('$N_{as}$')
        ax.set_ylabel('number')
#        ax.set_ylim(-100,ax.axis[3])
        fig.legend()
        #--- mean
#        print [m0,m1],list_of_zero_afs
#        pdb.set_trace()
        print 'n_ms=',len(df_n)#+len(df_zeros)
        print 'n_ms(0)=',len(df_zeros)
        N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
        err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
#        err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
    except:
        continue
    
#= [1.321e+03 6.500e+02 1.750e+02 4.800e+01 3.000e+01 1.300e+01 9.000e+00
 1.000e+00 1.000e+00 3.000e+00 3.000e+00 1.000e+00 1.000e+00 1.000e+00
 1.000e+00 1.000e+00 1.000e+00]
N= [  0.           0.89355469   1.78710938   2.97851562   3.87207031
   4.765625     5.95703125   6.85058594   7.74414062   8.93554688
   9.82910156  11.9140625   26.80664062  43.78417969  44.97558594
  68.80371094 304.70214844]
n_ms= 939
n_ms(0)= 1321
In [56]:
#--- choose bins for rho and lambda

import traceback

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)



#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()

#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )

#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']


#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)


#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
#ax.set_xlim(mc,swarm_lohi['magnitude'].max()+0.5)
m_min = min(m_list)-1.0e-6
m_max = max(m_list)+1.0e-6 #--- exclude the largest event
ax.set_xlim(m_min,np.ceil(m_max))
ax.set_ylim(0.5,1000)
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
           df_nn,
          'o',color='black',alpha=.1)
ax.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
         df_zeros,
          'x',color='gray')

plt.yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')

#--- choose m0<m<m1
dm = 0.75 #0.5
n = int( (6.0 - m_min ) / dm )
m = list(np.linspace(m_min,6.0,n))
#m=AdjustBins(m_list,m)
# m=[2.7999989999999997,
#  3.133332555555555,
#  3.466666111111111,
#  3.7999996666666664,
#  4.133333222222222,
#  4.466666777777778,
#  4.800000333333333,
#  5.800001]
m+=[6.41,7.1-0.01,7.1+0.01]
Mlist=m

m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
    try:
        m0 = m[ i ]
        m1 = m[ i + 1 ]
        df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
                                                         swarm_lohi,df_trig) #--- get parent with m0<m<m1
        df.reset_index( inplace = True, drop = True )
#        pdb.set_trace()
        #--- number of daughters for each mother
        df_n = df.groupby(by='Parent_id').count()['Event_id']
    
        #--- include daughters with zero children
        list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
        subset = swarm_lohi.loc[list_of_zero_afs]
        list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index    
        df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
        #--- mean
#        print [m0,m1],list_of_zero_afs
#        pdb.set_trace()
        
        N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
        err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
#        err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
    except:
#        traceback.print_exc()
        continue
    
#ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
            yerr=err_N,xerr=m_err,
            markerfacecolor='white',
            fmt='o', markersize=8,color='red',markeredgewidth=0.7,
            linewidth=.5,
            barsabove=None,capsize=5,capthick=1,elinewidth=1)


#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))

DrawFrame(ax, (0.15,0.08),(0.14,0.06),0.1,LOG_Y=True)

try:
    from scipy.optimize import curve_fit
    def func(x,c,alpha):
        return c*10**(alpha*x)
    popt, pconv = curve_fit(func,m_mean[m_mean>=mc],N_as[m_mean>=mc])
    
    mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
         swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
    mm = np.array(mm)
    mm=list(set(mm))
    mm.sort()
    mm=np.array(mm)
    ax.plot(mm,func(mm,*popt),'-.g')
    print popt
    print pconv**0.5
except:
    pass

plt.savefig('%s/prod.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')


fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
            yerr=err_N,xerr=m_err,
            markerfacecolor='white',
            fmt='o', markersize=8,color='red',markeredgewidth=0.7,
            linewidth=.5,
            barsabove=None,capsize=5,capthick=1,elinewidth=1)
ax.set_xlim(m_min,np.ceil(m_max))
ax.set_ylim(0.5,1000)
ax.set_yscale('log')
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
[1.54502169e-04 8.86558982e-01]
[[0.00013159        nan]
 [       nan 0.05237575]]
Out[56]:
[<matplotlib.axis.XTick at 0x7fe7495c0710>,
 <matplotlib.axis.XTick at 0x7fe74a30f2d0>,
 <matplotlib.axis.XTick at 0x7fe74975f750>,
 <matplotlib.axis.XTick at 0x7fe749619ed0>,
 <matplotlib.axis.XTick at 0x7fe74960e410>,
 <matplotlib.axis.XTick at 0x7fe74960e910>,
 <matplotlib.axis.XTick at 0x7fe74960ee10>]
In [57]:
m
Out[57]:
[1.999999,
 2.99999925,
 3.9999995000000004,
 4.99999975,
 6.0,
 6.41,
 7.09,
 7.109999999999999]
In [58]:
#--- density plots: rho(r)
#--- split based on the mainshock's magnitude

import numpy as np
import math

#def AdjustBins(data,bins):
#    hist, xedges=np.histogram(data, #--- histogram
#                              bins=bins)
#    tmp = [xedges[0]]
#    for i,j in zip(xedges[1:],hist): #--- expand bins 
#        if j < 2:
#            continue
#        tmp.append(i)
#    tmp.append(xedges[-1])
#    if hist[-1] < 2:
#        tmp.pop(-2)
#    return tmp 

#def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
#    df_parents = df.groupby(by='Parent_id').groups #--- parent events
#    ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
#    ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
#    df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
#    return df_m, np.mean( ds_m )

#def AddDist( df_trig, df_complete ):
#    x=df_complete.loc[ df_trig['Event_id'] ]['r(m)'] 
#    y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
#    df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
#    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )

def AddTime( df_trig, df_complete ):
    prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
    x=df_complete.loc[ df_trig['Event_id'] ]
    y=df_complete.loc[ df_trig['Parent_id'] ]
    x.reset_index(inplace=True)
    y.reset_index(inplace=True)
    df_trig['T'] = x['date']-y['date']
    df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)

def RemovePair( df, cs ):
    return df [ df[ 'R' ] <= df[ 'T' ] * cs ]

#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)

##--- parents' magnitudes
#df_trig = nij_trig.copy()
#m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#m_min = min(m_list)-1.0e-6
#m_max = max(m_list)+1.0e-6 #--- exclude the largest event
#
##--- choose m0<m<m1
#dm = 0.3
#n = int( ( m_max - m_min ) / dm )
#m = list(np.linspace(m_min,m_max,n))
#m=AdjustBins(m_list,m)
Mlist=m

#--- choose m0<m<m1
X={}
rho_dict={}
m_range={}
Err={}
m_mean={}
for i in xrange( len(m) - 1 ):
    try:
        m0 = m[ i ]
        m1 = m[ i + 1 ]
        
        df_trig, m_mean[ i ],junk = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
        df_trig.reset_index( inplace = True, drop = True )
        
        #--- add distance & time column
    #    AddDist( df_trig, df_complete )
    #    df_trig['R'] = rmat[df_trig['Parent_id'],df_trig['Event_id']]
    #    df_trig['T'] = tmat[df_trig['Parent_id'],df_trig['Event_id']]
    #    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
    #    AddTime( df_trig, df_complete )
        
        #--- remove pairs with r > cs.t
    #    cs = 3 * (24*3600) #--- km per day
    #    df_trig = RemovePair( df_trig, cs )
    
        #--- rho plots
        n_decades = int( math.ceil( log(df_trig['R'].max()/df_trig['R'].min(),10) ) )
        nbin_per_decade = 4
        nbins =  nbin_per_decade * n_decades
        
        bins = np.logspace(log(df_trig['R'].min(),10)-1.e-10,
                                                  log(df_trig['R'].max(),10)+1.e-10,nbins)
        bins = AdjustBins(df_trig['R'],bins)
    
        rho, xedges=np.histogram(df_trig['R'], 
                                 density=True,
                                 bins=bins)
        hist, xedges=np.histogram(df_trig['R'],
                                  bins=bins)
        x_bins, junk=np.histogram(df_trig['R'], #--- histogram
                                  bins=bins,
                                 weights = df_trig['R'])
    
        #--- rm empty bins
        x_bins = x_bins[ hist != 0 ]
        rho = rho[ hist != 0 ]
        hist = hist[ hist != 0 ]
        
        #--- poisson err 
        err = rho/np.sqrt(hist)
        x_bins /= hist
                
        #--- scattered plot
        fig=plt.figure(figsize=(8,4))
        ax=fig.add_subplot(1,2,1)
        ax.set_title('$%2.1f<m<%2.1f$'%(m0,m1))
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_xlabel('$t$')
        ax.set_ylabel('$r$')
        ax.set_xlim(1e-5,10)
        ax.set_ylim(1e-4,1e2)
        ax.scatter(df_trig['T'],df_trig['R'],alpha=0.1)
        
        #--- discretization effects
#         dt_min=df_trig['T'].min()
#         dt_max=df_trig['T'].max()
#         dr_min=df_trig['R'].min()
#         dr_max=df_trig['R'].max()
        
#         ax.plot([dt_min,dt_min],[dr_min,dr_max],
#                  '-',color='blue')
#         ax.plot([dt_max,dt_max],[dr_min,dr_max],
#                  '-',color='blue')
#         ax.plot([dt_min,dt_max],[dr_min,dr_min],
#                  '-',color='blue')
#         ax.plot([dt_min,dt_max],[dr_max,dr_max],
#                  '-',color='blue')
        
#         ax.set_ylim(dr_min/2,dr_max*2)
#         ax.plot([dt_min,dt_max],[dt_min*cs,dt_max*cs],
#                  '-',color='red')
        
        #--- rho
        ax=fig.add_subplot(1,2,2)
        ax.set_xlim(1e-4,1e2)
        ax.set_ylim(1e-5,1e1)
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_xlabel('$r$')
        ax.set_ylabel('$P(r)$')
        ax.errorbar(x_bins,rho, yerr=err)
            
        X[i] = x_bins[:]
        rho_dict[i]=rho[:]
        m_range[i]=np.array([m0,m1])
        Err[i]=err[:]
        plt.savefig('%s/rho_r.%s.png'%(DIR_OUTPT_figs,i),dpi=50,bbox_inches='tight')
    except:
        continue
In [59]:
#--- density plots: rho(r) 2nd


plt.figure(figsize=(15,5))

#--- plot rho
plt.subplot(1,3,1)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$r$')
plt.ylabel(r'$P(r)$')
plt.xlim(1e-3,1e2)
plt.ylim(1e-4,1e1)
for key in X:
#    if key < 3:
#        continue
    plt.errorbar(X[key],rho_dict[key],
                 yerr=Err[key],
                 fmt='-o', markersize=8
                )

#--- rescaled in x and y
plt.subplot(1,3,2)
alpha=0.5
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)
plt.xlim(1e-5,1e0)
plt.ylim(1e-2,1e3)

xdata = []
ydata = []
for key in X:
#    if key < 3:
#        continue
    xdata += list(X[key]/10**(alpha*np.mean(m_range[key])))
    ydata += list(rho_dict[key]*10**(alpha*np.mean(m_range[key])))
    
    plt.errorbar(X[key]/10**(alpha*m_mean[key]),
                 rho_dict[key]*10**(alpha*m_mean[key]),
                 yerr=Err[key]*10**(alpha*m_mean[key]),
                 label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
                 fmt='-o', markersize=8
                )
plt.legend(title='m',prop=fontP,loc='upper right',bbox_to_anchor=(0.4, .6))

#--- rescaled by r^-\nu
nu = 1.4 #2.9
plt.subplot(133)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$R$')
plt.ylabel(r'$P(R)\times R^{%s}$'%nu)
#plt.ylabel(r'$\lambda(t) \times t^p$')
plt.xlim(1e-4,1e0)
#plt.ylim(1e-3,1e5)
KEYS=X.keys()
KEYS.sort()
for key, index in zip(KEYS,xrange(100)):
#    if key < 3:
#        continue
    x_new = X[key]/10**(alpha*m_mean[key])
    y_new = rho_dict[key]*10**(alpha*m_mean[key])
    y2_new = Err[key]*10**(alpha*m_mean[key])
    plt.errorbar(x_new,
                 y_new*x_new**nu,
                 yerr=y2_new*x_new**nu, 
                 fmt='-o', markersize=8
                )
    
    
plt.tight_layout()
plt.savefig('%s/rho_r.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
In [60]:
#--- density plots: rho(r) (paper version)

import matplotlib as mpl
#mpl.rcParams.update(mpl.rcParamsDefault)
import matplotlib.ticker

#rc('text', usetex=True)
#font = {'size'   : 20}
#matplotlib.rc('font', **font)
#plt.style.use('seaborn')
#plt.rcParams.update({'lines.markeredgewidth': 1})

def returnCMtoInch(x,y):
    return (x/2.54,y/2.54)

fig, ax = plt.subplots(figsize=(4,4))
sigma=alpha
ax.axis([10/(100/1e-3),10,1e-3,100])
ax.loglog()

#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))

#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
#plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)


xdata = []
ydata = []
colors = ['black','red','green','blue','cyan','brown','violet','magenta','orange','yellow']
fillstyles=['white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None]
markers=['o','s','D','^','<','>','v']
assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
    xdata += list(X[key]/10**(sigma*np.mean(m_range[key])))
    ydata += list(rho_dict[key]*10**(sigma*np.mean(m_range[key])))
    
    ax.errorbar(X[key]/10**(sigma*m_mean[key]),
                 rho_dict[key]*10**(sigma*m_mean[key]),
                 yerr=Err[key]*10**(sigma*m_mean[key]),
                 label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
                 fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
                linewidth=.5,
                 barsabove=None,capsize=5,capthick=1,elinewidth=1
                )
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')

nu=-1.4
ax.plot([5e-3,.2],[10, 10* (.2/5e-3)**nu],'r-.',linewidth=1)

##--- rescaled by r^-\nu
#nu = 2.8
#plt.xscale('log')
#plt.yscale('log')
#plt.xlabel(r'$R$')
#plt.ylabel(r'$P(R)\times R^{%s}$'%nu)
##plt.ylabel(r'$\lambda(t) \times t^p$')

#KEYS=X.keys()
#KEYS.sort()
#for key, index in zip(KEYS,xrange(100)):
##    if key < 3:
##        continue
#    x_new = X[key]/10**(alpha*m_mean[key])
#    y_new = rho_dict[key]*10**(alpha*m_mean[key])
#    y2_new = Err[key]*10**(alpha*m_mean[key])
#    plt.errorbar(x_new,
#                 y_new*x_new**nu,
#                 yerr=y2_new*x_new**nu, 
#                 fmt='-o', markersize=8
#                )

    
    
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)


#plt.tight_layout()
plt.savefig('%s/rho_r.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')
In [61]:
#--- density plots: rho(r) (paper version 2nd)

#plt.style.use('seaborn')
#plt.rcParams.update({'lines.markeredgewidth': 1})

def returnCMtoInch(x,y):
    return (x/2.54,y/2.54)

fig, ax = plt.subplots(figsize=(4,4))
#---
ax.spines['right'].set_visible(False) #--- half open
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#---
ax.axis([40/(20/5e-4),40,5e-4,20])
ax.loglog()
#ax.set_xticks(np.logspace(-1,3,4))

#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))

#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())



ax.set_xlabel(r'$r(km)$')
ax.set_ylabel(r'$\rho(r)$')
#plt.xlim(1e-5,1e0)
#plt.ylim(1e-5,1e0)

#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
    xlabel = r'$%2.1f-%2.1f$'%(m_range[key][0],m_range[key][1])
    if key == 5:
        xlabel = r'$6.4$'
    if key == 7:
        xlabel = r'$7.1$'
    #    if key ==
    ax.errorbar(X[key],
                 rho_dict[key],
                 yerr=Err[key],
                 label=xlabel,
                 fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
                linewidth=.5,
                 barsabove=None,capsize=5,capthick=1,elinewidth=1
                )

    ax.legend(title=r'$m$',prop=fontP,loc='lower left',frameon=False,
              fontsize ='xx-small',markerscale=1,bbox_to_anchor=(1, 0))

    
nu=-1.4
ax.plot([1e-1,100],[50, 50* (100/1e-1)**nu],'r-.',linewidth=1)    
    
#ax.tick_params(axis='y',left=True, right=True,which='both')
#ax.tick_params(axis='x',bottom=True, top=True,which='both')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)


#plt.tight_layout()
plt.savefig('%s/rho_r.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
In [62]:
from scipy.optimize import curve_fit

def func(x,pref, xmean,sigma):
    return pref * np.exp(-((x-xmean)/sigma)**2)

fig=plt.figure(figsize=(8,4))
ax=fig.add_subplot(121)
ax.axis([2e-3,2e-3*1e1/1e-3,1e-3,1e1])
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$r$')
ax.set_ylabel('$P(r)$')


r_peak={}
mm_mean={}
#colors = ['black','red','green','blue','cyan','brown']
#assert len(colors) >= len(X)
p0 = (1.0,0.1,.01)
for key, col in zip(X,colors):
    try:
#        print np.mean(X[key])
        ax.errorbar(X[key],rho_dict[key], 
                     yerr=Err[key],
                    linestyle='None', color=col,
                     fmt='-o', markersize=8
                    )            
        popt, pcov = curve_fit(func, X[key], rho_dict[key],
                               sigma=Err[key], 
                               p0=p0) #, bounds=((0.0,0.0), (np.inf,10.0)))
        mm_mean[key]=m_mean[key]
        p0 = popt[:]
#        print p0
        r_peak[key] = popt[1]
    #    if key == 5:
        ax.plot(X[key], func(X[key], *popt),'-',color=col)#,label='fit: r=%3.2f' % (popt[1]))
        ax.legend()
    except:
        continue

ax=fig.add_subplot(122)
ax.set_yscale('log')
ax.set_xlabel('$m$')
ax.set_ylabel('$r$')
#ax.set_ylim(1,1e3)
ax.plot(mm_mean.values(),r_peak.values(),'o')
#
#
def func(x,c,alpha):
    return c*10**(alpha*x)
#
popt, pconv = curve_fit( func, mm_mean.values(),r_peak.values() )
#
ax.plot(mm_mean.values(),func( np.array(mm_mean.values()),*popt),'--r')
print popt
print pconv**0.5

fig.tight_layout()
[0.00095037 0.54468294]
[[0.02307428        nan]
 [       nan 1.51036767]]
In [63]:
#--- fitting

from scipy.optimize import curve_fit

def func(x,a,b):
    return a/x**b

xdata = np.array(xdata)
xxdata = xdata[xdata>2e-2]
ydata = np.array(ydata)
yydata = ydata[xdata>2e-2]

popt, pcov = curve_fit(func, xxdata, yydata, p0=(1.0,1.0), bounds=((0.0,0.0), (np.inf,10.0)))

plt.figure(figsize=(5,5))
plt.xscale('log')
plt.yscale('log')
#plt.xlim(1e-1,1e2)
#plt.ylim(1e-4,1e-1)
plt.plot(xdata, ydata,marker='o',linestyle='None')
plt.plot(xdata, func(xdata, *popt),'r-',label='fit: a=%3.2f' % popt[1])
plt.legend()         
#plt.plot(xdata, func(xdata, *popt), 'g--',
#         label='fit: a=%5.3f, b=%5.3f' % tuple(popt))
Out[63]:
<matplotlib.legend.Legend at 0x7fe748402a50>
In [64]:
#--- omori relation

import numpy as np
import math
from scipy.optimize import curve_fit

#def func(x,c,p,A):#,B):
#    return A/(1+x/c)**p

def func(x,c,p,A,c2,p2):#,B):
    return A/(1+(x/c)**p+(x/c2)**p2)


#def func(x,c,p):#,B):
#    return (p-1)/c/(1+x/c)**p



#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi.reset_index(inplace=True,drop=True)

popt={}
pcov={}
perr = {}

##--- parents' magnitudes
#m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#m_min = min(m_list)-1.0e-6
#m_max = max(m_list)+1.0e-6 #--- exclude the largest event
#
##--- choose m0<m<m1
#dm = 0.3
#n = int( ( m_max - m_min ) / dm )
#m = list(np.linspace(m_min,m_max,n))
#m=AdjustBins(m_list,m)
#Mlist=m



X={}
rho_dict={}
#m_range={}
Err={}
m_mean=np.array(np.zeros(len(m)-1))
Xedges={}
nrows=len(m)-1;ncols=2
fig=plt.figure(figsize=(ncols*4,nrows*4))


for i in xrange( len(m) - 1 ):
#    if i == (len(m) - 2):
#        pdb.set_trace()
    try:
        m0 = m[ i ]
        m1 = m[ i + 1 ]
        
        df_trig,m_mean[ i ],junk = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
        df_trig.reset_index( inplace = True, drop = True )
        print "# of mother events:%s"%len(df_trig.groupby(by="Parent_id"))
    
        #--- number of daughters for each mother
        df_n = df_trig.groupby(by='Parent_id').count()['Event_id']
        
        #--- include daughters with zero children
        list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
        subset = swarm_lohi.loc[list_of_zero_afs]
        list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index    
        df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
        #--- mean
        N_as = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
        print 'N_as=%s'%N_as
        #--- add distance & time column
    #    df_trig['R'] = rmat[df_trig['Parent_id'],df_trig['Event_id']]
    #    df_trig['T'] = tmat[df_trig['Parent_id'],df_trig['Event_id']]
    #    assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
        
        #--- remove pairs with r > cs.t
    #    cs = 3 * (24*3600) #--- km per day
    #    df_trig = RemovePair( df_trig, cs )
    
        #--- rho plots
        n_decades = int( math.ceil( log(df_trig['T'].max()/df_trig['T'].min(),10) ) )
        nbin_per_decade = 4
        nbins =  nbin_per_decade * n_decades
        bins = np.logspace(log(df_trig['T'].min(),10)-1.e-10,
                                                  log(df_trig['T'].max(),10)+1.e-10,nbins)
        bins = AdjustBins(df_trig['T'],bins)
    
        rho, xedges=np.histogram(df_trig['T'], #--- density
                                 density=True,
                                 bins=bins)
        rho *= N_as #--- int rho = N_as
        hist, xedges=np.histogram(df_trig['T'], #--- histogram
                                  bins=bins)
    
        x_bins, junk=np.histogram(df_trig['T'], #--- histogram
                                  bins=bins,
                                 weights = df_trig['T'])
    
        #--- rm empty bins
        x_bins = x_bins[ hist != 0 ]
        rho = rho[ hist != 0 ]
        hist = hist[ hist != 0 ]
        
        #--- poisson err 
        err = rho/np.sqrt(hist)
        x_bins /= hist
    
        #---- fit    
    #    fit = plfit.plfit( np.array(df_trig['T']) ) #--- clauset
    #    cval = fit[1] #--- completeness mag.
    #    pval = fit[0]
    #    [pval, gof] = plpva.plpva(np.array(df_trig['T']), cval, 'xmin', cval,'silent')
    #    print pval
        
        #--- 2nd fit
        try:
            popt[i], pcov[i] = curve_fit(func, x_bins, rho, 
                                   p0=(1.0e-2,1.0,1e2,1e0,2.0), 
                                   bounds=((1.0e-6,0.1,1e-5,1e-6,0.1), (np.inf,7,np.inf,np.inf,7)), #c,p,A,c2,p2
                                   sigma=err #rho**0.5 #err
        #                           absolute_sigma = True,
        #                          method='trf'#'dogbox'
                                  )
            perr[i] = np.sqrt(np.diag(pcov[i]))
        except:
            pass
        #--- scattered plot
        ax=fig.add_subplot(nrows,ncols,2*i+1)
        ax.set_title(r'$%2.1f<m<%2.1f$'%(m0,m1))
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_xlabel('$t$')
        ax.set_ylabel('$r$')
        #plt.xlim(1e-3,100)
        #plt.ylim(1e-3,1e2)
        ax.scatter(df_trig['T'],df_trig['R'],alpha=0.1)
        
        #--- discretization effects
        dt_min=df_trig['T'].min()
        dt_max=df_trig['T'].max()
        dr_min=df_trig['R'].min()
        dr_max=df_trig['R'].max()
        
        ax.plot([dt_min,dt_min],[dr_min,dr_max],
                 '-',color='blue')
        ax.plot([dt_max,dt_max],[dr_min,dr_max],
                 '-',color='blue')
        ax.plot([dt_min,dt_max],[dr_min,dr_min],
                 '-',color='blue')
        ax.plot([dt_min,dt_max],[dr_max,dr_max],
                 '-',color='blue')
        
        ax.set_ylim(dr_min/2,dr_max*2)
        ax.plot([dt_min,dt_max],[dt_min*cs,dt_max*cs],
                 '-',color='red')
        
        #--- rho
        ax=fig.add_subplot(nrows,ncols,2*i+2)
        ax.set_xlim(1e-4,1e2)
    #    ax.set_ylim(1e-3,1e3)
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_xlabel('$t$')
        ax.set_ylabel('$\lambda(t)$')
        ax.errorbar(x_bins,rho, 
                     yerr=err,
                    linestyle='None',
                     fmt='-o', markersize=8
                    )
        
        #--- plot the fit
    #    plt.plot(x_bins,
    #             (pval-1)/cval/(1+x_bins/cval)**pval,
    #             label='fit: p=%2.1f, c=%4.3f' % (pval,cval)
    #            )
        #--- 2nd fit
        try:
            ax.plot(x_bins, func(x_bins, *popt[i]),'r-',label='fit: p0=%3.2f, p1=%3.2f' % (popt[i][1],popt[i][4]))
        except:
            pass
        fig.tight_layout()
        ax.legend(prop=fontP)#bbox_to_anchor=(1.5, 0.5))
        
        X[i] = x_bins[:]
        Xedges[i]=xedges
        rho_dict[i]=rho[:]
    #    m_range[i]=np.array([m0,m1])
        Err[i]=err[:]
    except:
        continue
fig.savefig('%s/rho_t.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
# of mother events:594
N_as=0.6931927133269415
# of mother events:295
N_as=1.6023738872403561
# of mother events:44
N_as=5.340909090909091
# of mother events:4
N_as=21.25
# of mother events:1
N_as=69.0
# of mother events:1
N_as=305.0
In [65]:
KEYS=popt.keys()
KEYS.sort()
M_Mean=[]
pval=[]
cval = []
cerr = []
aval = []
aerr = []
mm=[]
pperr=[]
fig = plt.figure(figsize=(12,4))
ax=fig.add_subplot(131)
for ikey in KEYS: #--- exclude the last point
    pval.append( min(popt[ikey][1],popt[ikey][4]) )
    cval.append( min(popt[ikey][0],popt[ikey][3]) )
    aval.append( popt[ikey][2] )
    aerr.append(perr[ikey][2])
    index = min([popt[ikey][1],1],[popt[ikey][4],4])[1]
    pperr.append(perr[ikey][index])
    cerr.append(perr[ikey][index-1])
    M_Mean.append(m_mean[ikey])
ax.errorbar(M_Mean,pval,yerr=pperr,
            marker='o',color='black')
#ax.set_xlim(0,4)
ax.set_ylim(0,3)
ax.set_xlabel('$m$')
ax.set_ylabel('$p$')
ax.set_title('nth=%2.1e'%n_thresh)
#ax.plot([0,3],[1.29,1.29],'--r')
#ax.plot([0,3],[1.37,1.37],'--r')

#--- plot c 
def func(x,a,b):
    return a*10**(x*b)
poptt, pcovv = curve_fit(func, M_Mean,cval, sigma=cerr)
print poptt

ax2=fig.add_subplot(132)
ax2.set_ylim(1e-4,1e-1)
#ax2.set_xlim(0,4)
ax2.set_yscale('log')
ax2.set_ylabel('$c$')
ax2.set_xlabel('$m$')
ax2.errorbar(M_Mean,cval,yerr=cerr,
            marker='o',color='black')

try:
    ax2.plot(M_Mean,func(M_Mean,*poptt),'--r')
    
    #--- plot A
    poptt, pcovv = curve_fit(func, M_Mean,aval, sigma=aerr)
    ax3=fig.add_subplot(133)
    #ax3.set_ylim(1e2,1e4)
    #ax3.set_xlim(0,4)
    ax3.set_yscale('log')
    ax3.set_ylabel('$A$')
    ax3.set_xlabel('$m$')
    ax3.errorbar(M_Mean,aval,yerr=aerr,
                marker='o',color='black')
    ax3.plot(m_mean,func(M_Mean,*poptt),'--r')
    print poptt
except:
    pass

fig.tight_layout()
plt.savefig('%s/p.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
[2.30757697e-14 1.81870887e+00]
In [66]:
#--- omori relation


plt.figure(figsize=(16.5,5))


#--- unrescaled
plt.subplot(131)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$t$')
plt.ylabel('$\lambda(t)$')
#plt.xlim(1e-6,1e1)
#plt.ylim(1e0,1e7)
KEYS = X.keys()
KEYS.sort()
for key, index in zip(KEYS,xrange(100)):
    plt.errorbar(X[key],rho_dict[key],yerr=Err[key], 
#                label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),

                 fmt='-o', markersize=8
                )
#plt.legend(title='m',fontsize='small',loc='upper right',bbox_to_anchor=(0.4, .65),frameon=False)

#--- rescaled
alpha_c = 0.0 #0.5
alpha_p = 0.41 #--- productivity
alpha_y = alpha_c - alpha_p

plt.subplot(132)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$t/10^{%s\times m}$'%alpha_c)
plt.ylabel(r'$\lambda(t) \times 10^{(%s-%s)m}$'%(alpha_c,alpha_p))
#plt.ylabel(r'$\lambda(t) \times t^p$')
#plt.xlim(1e-6,1e0)
#plt.ylim(1e-3,1e5)
for key, index in zip(KEYS,xrange(100)):
#    exp=1.1 #pval[key]
    plt.errorbar(X[key]/10**(alpha_c*m_mean[key]),
                 rho_dict[key]*10**(alpha_y*m_mean[key]),
                 yerr=Err[key]*10**(alpha_y*m_mean[key]), 
#                 label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
                 label='%2.1f'%(m_mean[key]),
                 fmt='-o', markersize=8
                )
plt.legend(title='m',fontsize='small',loc='upper right',bbox_to_anchor=(0.4, .65),frameon=False)

#--- rescaled by t^-p
p = 1

plt.subplot(133)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$t/10^{%s\times m}$'%alpha_c)
plt.ylabel(r'$\lambda(t) \times 10^{(%s-%s)m} \times t^{%s}$'%(alpha_c,alpha_p,p))
#plt.ylabel(r'$\lambda(t) \times t^p$')
#plt.xlim(1e-6,1e0)
plt.ylim(1e-5,1e4)
for key, index in zip(KEYS,xrange(100)):
    x_new = X[key]/10**(alpha_c*m_mean[key])
    y_new = rho_dict[key]*10**(alpha_y*m_mean[key])
    y2_new = Err[key]*10**(alpha_y*m_mean[key])
    plt.errorbar(x_new,
                 y_new*x_new**p,
                 yerr=y2_new*x_new**p, 
                 fmt='-o', markersize=8
                )

plt.tight_layout()
plt.savefig('%s/rho_t_mult.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
In [67]:
#--- omori (paper version)


#--- rescaled

fig, ax = plt.subplots(figsize=(4,4))
ax.axis([4/(1000/3e-4),4,3e-4,1000])
ax.loglog()

##--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))

#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())


#xdata = []
#ydata = []
#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
#    xdata += list(X[key]/10**(alpha*np.mean(m_range[key])))
#    ydata += list(rho_dict[key]*10**(alpha*np.mean(m_range[key])))
    
    ax.errorbar(X[key]/10**(alpha_c*m_mean[key]),
                 rho_dict[key]*10**(alpha_y*m_mean[key]),
                 yerr=Err[key]*10**(alpha_y*m_mean[key]), 
                 label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
                 fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
                linewidth=.5,
                 barsabove=None,capsize=5,capthick=1,elinewidth=1
                )
        
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')

p=-1
ax.plot([1e-4,2],[5, 5* (2/1e-4)**p],'r-.',linewidth=1)

DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)

#plt.tight_layout()
plt.savefig('%s/rho_t.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')

$\lambda(t)$ vs. $t$

In [68]:
#--- density plots: rho(r) (paper version 2nd)


def returnCMtoInch(x,y):
    return (x/2.54,y/2.54)

fig, ax = plt.subplots(figsize=(4,4))
#---
ax.spines['right'].set_visible(False) #--- half open
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#---
ax.axis([10/(6e3/2e-3),10,2e-3,6e3])
ax.loglog()

#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))

#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

ax.set_xlabel(r'$t(day)$')
ax.set_ylabel(r'$\lambda(t)$')
#plt.xlim(1e-5,1e0)
#plt.ylim(1e-5,1e0)

#colors = ['black','red','green','blue','cyan']
#fillstyles=['white',None,'white',None,'white',None]
#markers=['o','s','D','^','<']
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
    xlabel = r'$%2.1f-%2.1f$'%(m_range[key][0],m_range[key][1])
    if key == 5:
        xlabel = r'$6.4$'
    if key == 7:
        xlabel = r'$7.1$'
    ax.errorbar(X[key],
                 rho_dict[key],
                 yerr=Err[key],
                 label=xlabel,
                 fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
                linewidth=.5,
                 barsabove=None,capsize=5,capthick=1,elinewidth=1
                )

#    ax.legend(title=r'$m$',prop=fontP,loc='lower left',frameon=False, 
#              bbox_to_anchor=(1, 0))
#ax.tick_params(axis='y',left=True, right=True,which='both')
#ax.tick_params(axis='x',bottom=True, top=True,which='both')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)


#plt.tight_layout()
plt.savefig('%s/rho_t.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)

rescaled $\lambda(t)$ vs. $t$

In [69]:
#--- density plots: rho(r) (paper version 2nd)


def returnCMtoInch(x,y):
    return (x/2.54,y/2.54)

fig, ax = plt.subplots(figsize=(4,4))
#---
ax.spines['right'].set_visible(False) #--- half open
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#---
ax.axis([10/(6e3/2e-3),10,2e-3,6e3])
ax.loglog()

#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))

#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())

#--- rescaled
alpha_c = 0.0 #0.5
alpha_p = 0.6 #--- productivity
alpha_y = alpha_c - alpha_p
p = 1

plt.xlabel(r'$t(day)$')
plt.ylabel(r'$\lambda(t) \times t^{%s}$'%(p))

#plt.xlim(1e-5,1e0)
#plt.ylim(1e-5,1e0)

#colors = ['black','red','green','blue','cyan']
#fillstyles=['white',None,'white',None,'white',None]
#markers=['o','s','D','^','<']
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
    xlabel = r'$%2.1f-%2.1f$'%(m_range[key][0],m_range[key][1])
    if key == 5:
        xlabel = r'$6.4$'
    if key == 7:
        xlabel = r'$7.1$'
    ax.errorbar(X[key],
                 rho_dict[key]*X[key]**p,
                 yerr=Err[key]*X[key]**p, 
                 label=xlabel,
                 fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
                linewidth=.5,
                 barsabove=None,capsize=5,capthick=1,elinewidth=1
                )

    ax.legend(title=r'$m$',prop=fontP,loc='lower left',frameon=False, 
              bbox_to_anchor=(1, 0))
#ax.tick_params(axis='y',left=True, right=True,which='both')
#ax.tick_params(axis='x',bottom=True, top=True,which='both')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)


#plt.tight_layout()
plt.savefig('%s/rho_t_rescaled.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
In [70]:
#--- construct trees

from math import *

class tree():
    leafDepth = 0
    n_leaves = 0
    g_netx = None
    leaf_generation = {} #--- key: leaf id;  val: generation
    nsize = 1 #--- cluster size
    rtm_dic = {}
    #---
    def __init__( self, data, coord, coord_parent, n_generation ):
        self.data = data #--- id
        self.coord = coord
        self.coord_parent = coord_parent
        self.n_daughters = 0
        self.daughter_obj = {}
        self.n_generation = n_generation
    #---
    def insert( self, data, parent_id, coord, coord_parent, n_generation ):
        if parent_id == self.data: #--- parent exists
            self.daughter_obj[ self.n_daughters ] = tree( data, coord, coord_parent, n_generation )
            self.n_daughters += 1
            return
        else:
            for keys in self.daughter_obj:
                self.daughter_obj[ keys ].insert( data, parent_id, coord, coord_parent, n_generation )	
    #---
    def sprint( self, df ):
        for keys in self.daughter_obj:
            tree.g_netx.add_node(self.daughter_obj[ keys ].data)
            tree.g_netx.add_edge( self.data, self.daughter_obj[ keys ].data)
            self.daughter_obj[ keys ].sprint( df )
    #---
    def getLeaves( self ):
        if len( self.daughter_obj ) == 0: #--- it's a leaf
            tree.leaf_generation[ self.data ] = self.n_generation
        else:
            for keys in self.daughter_obj: #--- search in daughters
                self.daughter_obj[ keys ].getLeaves()
    #---
    def getSize( self ):
        for keys in self.daughter_obj: #--- search in daughters
            tree.nsize += 1
            self.daughter_obj[ keys ].getSize()
    
    def get( self, title ):
        pass


    def get_coords( self ): #--- rtm_dic['node_id'] = [r,t,m]
        for keys in self.daughter_obj:
            child_object = self.daughter_obj[keys]
            cords = child_object.coord
            id_node = child_object.data
            tree.rtm_dic[ id_node ] = cords
            child_object.get_coords()
            
    def plot(self,ax):
#        print 'id=%s'%(self.data), 'generation=%s'%(self.n_generation)
        for keys in self.daughter_obj:
            child_object = self.daughter_obj[keys]
            cords = child_object.coord
            tree.pairs.append(child_object.data) #[self.coord, child_object.coord])
#            ax.annotate("", tuple(cords[1][:2]), xytext=tuple(self.coord[1][:2]),
#                         textcoords='data', annotation_clip=False,
#                        arrowprops=dict(arrowstyle="-|>,head_width=.1,head_length=.4",color="b",linewidth="0.3"))  
            child_object.plot(ax)

    def nodeExist(self, event_id):
        if self.data == event_id:
#            print 'yes'
            return True
        
        for keys in self.daughter_obj:
            child_object = self.daughter_obj[keys]
            if child_object.nodeExist( event_id ):
                return True
        return False

pref = 1.0 / (24.0*60*60)
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
                             'magnitude', 
                             ( mc, sys.maxint ) ) 
swarm_lohi = GetCartesian( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)

NoTriggerList = set(nij_trig['Parent_id'])-set(nij_trig['Event_id']) #--- list of roots
    
#---
#--- include parent and entire generation
#---
CLUSTER={}
clusterID = 0
stree = {}
magScale = 1.0
for rootID in NoTriggerList: #--- roots: events with no parents
    tree_level = 0
    mainShockList = [ rootID ]
    CLUSTER[clusterID]=[mainShockList[0]]
    t_root = swarm_lohi.loc[ rootID ]['date']
    r_root = swarm_lohi.loc[ rootID ][['x(km)','y(km)','z(km)']].tolist()
    m_root = swarm_lohi.loc[ rootID ]['magnitude'] 
    stree[ clusterID ] = tree( rootID, [ t_root, r_root, m_root ],
                                           [ t_root, r_root, m_root ], 
                                           tree_level )
    while len(mainShockList): #--- from root down to leaves
        newMainShockList = []
        for mainShockID, counter in zip( mainShockList, xrange( sys.maxint ) ):
            try:
                newMainShockList += d_trig[ mainShockID ] #--- aftershocks become new
            except:
#                traceback.print_exc()
                continue
            for aftershockID in d_trig[ mainShockID ]:
                CLUSTER[ clusterID ].append(aftershockID)
#                dt = (swarm_lohi.loc[ aftershockID ]['date']-swarm_lohi.loc[ rootID ]['date']).total_seconds() * pref 
                t_afs = swarm_lohi.loc[ aftershockID ]['date']
#                -swarm_lohi.loc[ rootID ]['date']).total_seconds() * pref 
                r_afs = swarm_lohi.loc[aftershockID][['x(km)','y(km)','z(km)']].tolist()
                m_afs = swarm_lohi.loc[aftershockID]['magnitude']
                t_parent = swarm_lohi.loc[ mainShockID ]['date']
                r_parent = swarm_lohi.loc[mainShockID][['x(km)','y(km)','z(km)']].tolist()
                m_parent = swarm_lohi.loc[mainShockID]['magnitude']
                stree[ clusterID ].insert( aftershockID, mainShockID, 
                                          [t_afs, r_afs, m_afs],
                                          [t_parent, r_parent, m_parent],
                                          tree_level + 1 )
        mainShockList = [ i for i in newMainShockList ]
        tree_level += 1
    assert len(mainShockList) == 0 #--- leaf
    clusterID+=1
    
sl=[]
for keys in stree:
    sl.append([stree[keys].n_daughters,stree[keys].data,keys])
sl.sort(reverse=True)
sl[:3]
Out[70]:
[[69, 3, 2], [27, 498, 293], [10, 1869, 205]]
In [71]:
x=(-97,-96.5)
y=(35.8,36)
swarm_lohi.sort_values(by='magnitude',ascending=False).head() #[(swarm_lohi['longitude']>x[0]) & (swarm_lohi['longitude']<x[1]) & (swarm_lohi['latitude']>y[0]) & (swarm_lohi['latitude']<y[1])].head(n=10)
Out[71]:
date latitude longitude depth magnitude ID x(km) y(km) z(km)
666 2019-07-06 03:19:52.860 35.771818 -117.594157 3.241 7.10 70611993 24.872893 32.932561 3.241
3 2019-07-04 17:33:48.760 35.708350 -117.498625 15.054 6.40 70463229 35.507467 25.890521 15.054
723 2019-07-06 03:47:53.460 35.904932 -117.745199 6.039 5.50 70613673 8.058975 47.702359 6.039
797 2019-07-06 04:18:55.800 35.904789 -117.683976 13.805 5.44 70615536 14.874288 47.686492 13.805
669 2019-07-06 03:23:50.370 35.790344 -117.598494 12.817 5.37 70612231 24.390100 34.988113 12.817
In [72]:
#--- a large event belongs to what cluster

event_id = 5 #1509
FOUND = False
for item in sl:
    cls_id = item[2]
    root_id = item[1]    
    if stree[cls_id].nodeExist(event_id):
        FOUND = True
        break
if FOUND:
    print 'root_id=%s'%root_id
else:
    print 'it is a background event'
root_id=3
In [73]:
#--- leaf vs. depth

import networkx as nx
#import treeDataStruct as tds
from itertools import count
import random as rnd

exp = 4
for item, index in zip(sl[:6],xrange(sys.maxint)):
    sid = item[2]
    root_id = item[1]
 #   print 'root_id=',root_id
    tree.g_netx = None
    tree.g_netx = nx.DiGraph(directed=True)
    tree.g_netx.add_node(stree[ sid ].data,color='b')
    stree[ sid ].sprint(swarm_lohi )
    
    # read edgelist from grid.edgelist
    nx.write_edgelist(tree.g_netx, path="grid.edgelist", delimiter=":")
    H = nx.read_edgelist(path="grid.edgelist", delimiter=":") #,create_using=nx.DiGraph())
#    nx.draw(H,with_labels=True,arrows = True)#pos=nx.spring_layout(H))
#    nx.draw_networkx(tree.g_netx, arrows=True, with_labels=True,node_color='b') #, **kwds)

    pos = nx.spring_layout(H)
    nodes=H.nodes

    colors = ['red' for n in nodes]
    sizes=[exp**swarm_lohi.loc[ int(n) ]['magnitude'] for n in nodes]

    for n,indx in zip(nodes,xrange(sys.maxint)):
        if int(n) == root_id:
            colors[indx] = 'gray'
    ec = nx.draw_networkx_edges(H, pos, alpha=0.2,arrows = True)
    nc = nx.draw_networkx_nodes(H, pos,node_size=sizes, node_color=colors,cmap=plt.cm.seismic) #, nodelist=nodes, , 
 #                           with_labels=False, node_size=100, cmap=plt.cm.jet)

    plt.axis('off')
    plt.savefig('%s/topology.%s.png'%(DIR_OUTPT_figs,index),dpi=75)#,bbox_inches='tight')
    plt.show()

size vs. depth

In [74]:
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings

rc('text', usetex=True)
font = {'size'   : 26}
matplotlib.rc('font', **font)

fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xscale('linear')
ax.set_yscale('log')
#ax.set_xlabel('mean depth')
#ax.set_ylabel('size')
ax.axis([0,7.0,1,1e2])
ax.set_yscale('log')
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(ax.axis()[0],ax.axis()[1],int(nbin)+1))

#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))

MD=[]
S=[]
M=[]
root=[]
for item in sl:
    sid = item[2]
    root_id = item[1]
    root.append(root_id)

    tree.leaf_generation={} #--- initialize
#    stree[ sid ].initialize()
    stree[ sid ].getLeaves() #--- call function
    d = stree[ sid ].leaf_generation
#    print 'd[%s]=%s'%(sid,d)
    mean_depth = np.mean( d.values() ) #np.mean( d.values() ) #1.0 * tds.tree.leafDepth / tds.tree.n_leaves
    #--- size
    tree.nsize = 1 #--- initialize
    stree[ sid ].getSize() #--- call function
    MD.append(mean_depth)
    S.append(stree[ sid ].nsize)
    M.append(swarm_lohi.loc[stree[ sid ].data]['magnitude'])

exp = 5*.7
S=np.array(S)
MD = np.array(MD)
M=np.array(M)
root=np.array(root)
scatter=ax.scatter(MD, S,
            marker='o',color='gray',edgecolors='none',alpha=0.5,
            s=exp**np.array(M))

#--- color specific clusters
cond = np.all([S>=5,MD>0.0],axis=0)
scatter=ax.scatter(MD[cond], S[cond],
            marker='o',color='red',edgecolors='none',#alpha=0.5,
            s=exp**np.array(M[cond]))
#print root[10<=S], S[10<=S] 
print 'roots=',root[cond]
print 'size',S[cond]
print 'depth',MD[cond]

# #--- legend
# M=np.array(M)
# m_min=M.min()
# m_max=M.max()
# m_mean=(m_min+m_max)/2 #M.mean()
# l1 = plt.scatter([-2],[-2], s=exp**m_min, c='gray',edgecolors='none')
# l2 = plt.scatter([-3],[-3], s=exp**m_mean,c='gray', edgecolors='none')
# l3 = plt.scatter([-4],[-4], s=exp**m_max,c='gray', edgecolors='none')
# labels = ["%2.1f"%m_min, "%2.1f"%m_mean, "%2.1f"%m_max]
# leg = plt.legend([l1, l2, l3], labels, ncol=3,fontsize=16, frameon=False,
#                 handlelength=0, borderpad = 0, labelspacing=0,columnspacing=2,prop=None,
#                 handletextpad=1, scatterpoints = 1,loc='upper right', bbox_to_anchor=(.92, 0.95))
#--- frame
DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.1,LOG_Y=True)

#--- save
plt.savefig('%s/topology.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')

mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rc('text', usetex=True)
font = {'size'   : 26}
matplotlib.rc('font', **font)
roots= [   3  498 1869   71 3610 1898 1670 2494 2425 2213 1389  357  828  640
  351  335  305 3635 3503 2787 2108 1930 1665 1643 1343 1286  848  768
  663  391  129 3845 3543 1492]
size [ 254   48   17    9    6    9   13   21    9    6    6    6    8    5
    8    9    6    5    5    5    5    5    9    7    5    5    6   10
 1071    5   11    7   18    5]
depth [3.11643836 1.82857143 1.6        1.33333333 1.         1.6
 3.         3.5        2.         1.25       1.25       1.25
 2.         1.33333333 2.33333333 2.4        1.66666667 1.66666667
 2.         1.66666667 2.         2.         2.75       2.
 1.66666667 2.         1.75       2.6        8.62865497 2.
 2.66666667 2.         3.6        2.5       ]

trig. cascades

In [75]:
#--- leaf vs. depth

import networkx as nx
#import treeDataStruct as tds
from itertools import count
import random as rnd



#dd={'ROOT':407,
#    'exponent':4.0*2.7}
#dd={'ROOT':325,
#      'exponent':4.0*2.7}
dd={'ROOT':1,#1505,
     'exponent':4.0*.8}

#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']

for item in sl:
    cls_id = item[2]
    root_id = item[1]
    if root_id != ROOT:
        continue
    tree.pairs = []
    stree[cls_id].plot(ax)
    tree.pairs = [root_id]+tree.pairs
    
for item in sl:
    sid = item[2]
    root_id = item[1]
    if root_id != ROOT:
        continue
    print 'root_id=',root_id, 'n=%s'%(item[0])
    tree.g_netx = None
    tree.g_netx = nx.DiGraph(directed=True)
    tree.g_netx.add_node(stree[ sid ].data,color='b')
    stree[ sid ].sprint(swarm_lohi )
    
    # read edgelist from grid.edgelist
    nx.write_edgelist(tree.g_netx, path="grid.edgelist", delimiter=":")
    H = nx.read_edgelist(path="grid.edgelist", delimiter=":") #,create_using=nx.DiGraph())
#    nx.draw(H,with_labels=True,arrows = True)#pos=nx.spring_layout(H))
#    nx.draw_networkx(tree.g_netx, arrows=True, with_labels=True,node_color='b') #, **kwds)

    pos = nx.spring_layout(H)
    nodes=H.nodes

    colors = ['red' for n in nodes]
    
    sizes=[exponent**swarm_lohi.loc[ int(n) ]['magnitude'] for n in nodes]
    node_shape = ['o' for n in nodes]
    for n,indx in zip(nodes,xrange(sys.maxint)):
        if int(n) == root_id:
            colors[indx] = 'yellow'
#            node_shape[indx]='*'

    fig=plt.figure(figsize=(8,8))
    ax=fig.add_subplot(111)
    ec = nx.draw_networkx_edges(H, pos, alpha=0.2,arrows = True, ax= ax,edge_color='b',width=3)
    nc = nx.draw_networkx_nodes(H, pos,node_size=sizes, node_color=colors,cmap=plt.cm.seismic, ax=ax,
                               edgecolors='black',node_shape='o') #, nodelist=nodes, , 
 #                           with_labels=False, node_size=100, cmap=plt.cm.jet)

    plt.axis('on')
    ax.tick_params(labelbottom=False,labelleft=False,length=0.0)
    fig.show()
#    DrawFrame(ax, (0.0,0.0),(0.0,0.0),.1)
#    DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
    DrawFrame(ax, (0.16,0.1),(0.09,0.03),0.1)


#plt.tight_layout()
    fig.savefig('%s/tree.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
    
    

timeseries associated with the trig. cascade

In [76]:
#--- plot clusters


def Inside(t,(tmin,tmax)):
    if tmin<= t<tmax:
        return True
    return False    

def GetTrigList( nij_trig ):
    d_trig = {}
    for items in nij_trig.itertuples():
        triggerID = items.Parent_id
        d_trig.setdefault( triggerID, [] ).append( items.Event_id )
    return d_trig

def PlotArrowsTime(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
    tlist=[]
    mlist=[]
    xm = []
    ym = []
    for triggerID in d_trig: #--- mother events
        
        t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
        m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
        if not ( Inside(t0,(tmin,tmax) ) and #--- within the given interval
                 Inside(swarm_lohi['x(km)'].iloc[triggerID],(xmin,xmax)) and 
                 Inside(swarm_lohi['y(km)'].iloc[triggerID],(ymin,ymax)) ):
            continue

#        x0 = md.date2num( t0 ) #--- convert to time object
#        y0 = m0 #--- magnitude
        
#        tlist.append( t0 )
#        mlist.append( m0 )

        for daughter_id in d_trig[triggerID]: #--- children
            
            t1 = df_complete['date'].iloc[ daughter_id ] #--- time
            m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
            
            
            if not ( Inside(t1,(tmin,tmax) ) and #--- within the given interval
                     Inside(swarm_lohi['x(km)'].iloc[daughter_id],(xmin,xmax)) and 
                     Inside(swarm_lohi['y(km)'].iloc[daughter_id],(ymin,ymax)) ):
                continue
                
#            x1 = md.date2num( t1 ) #--- convert to time object
#            y1 = m1 #--- magnitude

#            xw = x1 - x0 #--- t_child-t_parent
#            yw = y1 - y0 #--- m_child - m_parent
            if daughter_id in tree.pairs:
                tlist.append( t1 )
                mlist.append( m1 )
                ax.annotate("", (t1,m1), xytext=(t0, m0), zorder=1,
                        textcoords='data',clip_on=True,
                        arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
                                        shrinkA=shrinkA,
                                        color="b",linewidth="0.3")) 
    
        #--- plot mother
        if triggerID == tree.pairs[0]: #in roott:
            xm.append(t0)
            ym.append(m0)
            
            #--- plot circles
    df=pd.DataFrame({'date':tlist,'mag':mlist})
    ax.scatter(df['date'], df['mag'],
                s=exponent**(df['mag']),
                alpha=1,zorder=2,
                facecolors='red',color='black')
    #plt.savefig('timeSeries.png')
    for x0,y0 in zip(xm,ym):
#        print x0,y0
        ax.scatter(x0,y0,
                    s=exponent**y0,
                    marker='*', zorder=3,
                    facecolors='yellow',color='black',
                    alpha=1.0)
        
#--------------------------------------
#----- key: event value: aftershock id
#--------------------------------------
d_trig = GetTrigList( nij_trig )

# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# swarm_lohi.reset_index(inplace=True,drop=True)

# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )



# dd={'ROOT':407,
#     'exponent':4.5*3.2,
#     'dt_before':0.005,'dt_after':.1,
#     'head_width':.2,'head_length':.4,
#       'shrinkA':0}
# dd={'ROOT':325,
#     'exponent':4.5*3.2,
#     'dt_before':0.005,'dt_after':.1,
#     'head_width':.2,'head_length':.4,
#     'shrinkA':0} 
dd={'ROOT':1,#1505,
    'exponent':4.5*.8,
    'dt_before':0.01,'dt_after':0.25*12,
    'head_width':.2,'head_length':.4,
     'shrinkA':0}

#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']
dt_before = dd['dt_before']
dt_after=dd['dt_after']
head_width = dd['head_width']
head_length=dd['head_length']
shrinkA=dd['shrinkA']
swarm_lohi = GetCartesian( swarm_lohi )
dff = swarm_lohi.sort_values('magnitude',ascending=False)
roott = [i[1] for i in sl]
#for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
for item in dff.itertuples(): #--- plot clusters associated with large events
    if item.Index != ROOT:
        continue
    
    (xc,yc)=(item.longitude,item.latitude) #(-96.807159,35.985256) #(-98.739380,36.467762) #(-99.0,36.5)
    tc = item.date #pd.to_datetime("2016-09-03 12:02:44.520")
    dx = 1000

    fig = plt.figure(figsize=(24,8),frameon=False) #, dpi=150)
    ax = fig.add_subplot(111)
    
    #--- plot all events
    plt.ylim(mc,swarm_lohi['magnitude'].max())
    plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
                s=exponent**swarm_lohi['magnitude'],
            alpha=0.1)

    tt = tc #swarm_lohi['date'].loc[key_max]
    xx={0:xc,1:yc}

    PlotArrowsTime( swarm_lohi, d_trig,
                       ( tt+datetime.timedelta(days=-dt_before),tt+datetime.timedelta(days=+dt_after)), 
           (xx[0]-dx,xx[0]+dx),
           (xx[1]-dx,xx[1]+dx) )
#    ax.set_xlim(tt+datetime.timedelta(days=-dt_before),swarm_lohi.iloc[tree.pairs]['date'].max())
    ax.set_xlim(tt+datetime.timedelta(days=-dt_before),tt+datetime.timedelta(days=+dt_after))
    ax.set_ylim(mc-0.1,7.5)
    
    ax.tick_params(axis='y',labelsize=32,right=True)
    ax.tick_params(axis='x',rotation=20,labelsize=28,top=True)
    ax.xaxis.set_major_formatter(dates.DateFormatter('%d %b %Y \n %H:%M'))
    #--- frame
    DrawFrame(ax, (0.04,0.04),(0.22,0.04),0.1)
    fig.savefig('%s/timeSeries.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)

TypeErrorTraceback (most recent call last)
<ipython-input-76-da80ca029a35> in <module>()
    143     plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
    144                 s=exponent**swarm_lohi['magnitude'],
--> 145             alpha=0.1)
    146 
    147     tt = tc #swarm_lohi['date'].loc[key_max]

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/pyplot.pyc in scatter(x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, hold, data, **kwargs)
   3376                          vmin=vmin, vmax=vmax, alpha=alpha,
   3377                          linewidths=linewidths, verts=verts,
-> 3378                          edgecolors=edgecolors, data=data, **kwargs)
   3379     finally:
   3380         ax._hold = washold

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1715                     warnings.warn(msg % (label_namer, func.__name__),
   1716                                   RuntimeWarning, stacklevel=2)
-> 1717             return func(ax, *args, **kwargs)
   1718         pre_doc = inner.__doc__
   1719         if pre_doc is None:

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, **kwargs)
   4021             linewidths = rcParams['lines.linewidth']
   4022 
-> 4023         offsets = np.column_stack([x, y])
   4024 
   4025         collection = mcoll.PathCollection(

/global/software/anaconda/anaconda-2.7-5.1.0/lib/python2.7/site-packages/numpy/lib/shape_base.pyc in column_stack(tup)
    367             arr = array(arr, copy=False, subok=True, ndmin=2).T
    368         arrays.append(arr)
--> 369     return _nx.concatenate(arrays, 1)
    370 
    371 def dstack(tup):

TypeError: invalid type promotion

spatial map associated with the trig. cascade

In [77]:
#--- plot clusters


def Inside(t,(tmin,tmax)):
    if tmin<= t<tmax:
        return True
    return False    

def GetTrigList( nij_trig ):
    d_trig = {}
    for items in nij_trig.itertuples():
        triggerID = items.Parent_id
        d_trig.setdefault( triggerID, [] ).append( items.Event_id )
    return d_trig

    
def PlotArrowsSpace(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
    xlist=[]
    ylist=[]
    mlist=[]
    xm=[]
    ym=[]
    m_mother=[]
    for triggerID in d_trig: #--- mother events
        
        t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
        ( x0, y0 ) = ( swarm_lohi['longitude'].iloc[triggerID], 
                       swarm_lohi['latitude'].iloc[triggerID] )
        m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
        if not ( #Inside(t0,(tmin,tmax) ) and #--- within the given interval
                  Inside(x0,(xmin,xmax)) and 
                  Inside(y0,(ymin,ymax)) ):
             continue

        for daughter_id in d_trig[triggerID]: #--- children
            t1 = df_complete['date'].iloc[ daughter_id ] #--- time
            ( x1, y1 ) = ( swarm_lohi['longitude'].iloc[ daughter_id ], 
                           swarm_lohi['latitude'].iloc[ daughter_id ] )
            m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
                        
            if not ( #Inside(t1,(tmin,tmax) ) and #--- within the given interval
                      Inside(x1,(xmin,xmax)) and 
                      Inside(y1,(ymin,ymax)) ):
                 continue

            if daughter_id in tree.pairs:
                xlist.append( x1 )
                ylist.append( y1 )
                mlist.append( m1 )
                ax.annotate("", (x1,y1), xytext=(x0, y0), zorder=1,
                         textcoords='data', annotation_clip=False, 
#                         arrowprops=dict(width=.2,headlength=.4,shrink=0.5,color='b'))
                            # backgroundcolor="w",
                        arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
                                        shrinkA=shrinkA,
                                        color="b",linewidth="0.3")) 
 
        #--- plot mother
        if triggerID == tree.pairs[0]: # in roott:
            xm.append(x0)
            ym.append(y0)
            m_mother.append(m0)
        

    #--- plot circles
    df=pd.DataFrame({'x(km)':xlist,'y(km)':ylist,'m':mlist})
    ax.scatter(df['x(km)'], df['y(km)'],
                s=exponent**(df['m']),
                alpha=1, zorder=2,
                facecolors='red',color='black')
    
    for x0,y0,m0 in zip(xm,ym,m_mother):
        plt.scatter(x0,y0,
                    s=exponent**m0,
                    marker='*',zorder=3,
                    facecolors='yellow',color='black',
                    alpha=1.0)
    
# #--------------------------------------
# #----- key: event value: aftershock id
# #--------------------------------------
# d_trig = GetTrigList( nij_trig )

# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# swarm_lohi.reset_index(inplace=True,drop=True)

# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )


# dd={'ROOT':407,
#     'exponent':4.0*2.7,
#     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#     swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}
# dd={'ROOT':325,
#     'exponent':4.0*2.7,
#      'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}
dd={'ROOT':1505,
    'exponent':4.0*0.8,
     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
    'head_width':.15,'head_length':.3,
     'shrinkA':0}

#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']
dx = dd['dx']
head_width = dd['head_width']
head_length=dd['head_length']
shrinkA=dd['shrinkA']
#--- plot arrows
#for ii in xrange(1):
swarm_lohi = GetCartesian( swarm_lohi )
dff = swarm_lohi.sort_values('magnitude',ascending=False)
roott = [i[1] for i in sl]
#for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
for item in dff.itertuples(): #--- plot clusters associated with large events
    if item.Index != ROOT:
        continue
    
    (xc,yc)=(item.longitude,item.latitude) #(-96.807159,35.985256) #(-98.739380,36.467762) #(-99.0,36.5)
    tc = item.date #pd.to_datetime("2016-09-03 12:02:44.520")
    dt = 1000   #--- setup

    
    fig = plt.figure(figsize=(8,8)) #,frameon=False) #, dpi=300)
    ax = fig.add_subplot(111)
    
    #--- plot all events
    ax.scatter(swarm_lohi['longitude'], swarm_lohi['latitude'], # 'x(km)'],swarm_lohi['y(km)'],
                s=exponent**swarm_lohi['magnitude'],
                alpha=0.1)

    tt = tc #swarm_lohi['date'].loc[key_max]
    xx={0:xc,1:yc}
    PlotArrowsSpace( swarm_lohi, d_trig,
    ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)), 
           (xx[0]-dx,xx[0]+dx),
           (xx[1]-dx,xx[1]+dx) )

    
    ax.set_xlim(xx[0]-dx,xx[0]+dx)
    ax.set_ylim(xx[1]-dx,xx[1]+dx)
    
    
    #--- set tick labels
    xlist=np.linspace(xx[0]-dx,xx[0]+dx,4,endpoint=True)
    ax.set_xticks(xlist)
    xlist=[r'$%5.3f^{\circ}$'%i for i in xlist]
    #
    ylist=np.linspace(xx[1]-dx,xx[1]+dx,4,endpoint=True)
    ax.set_yticks(ylist)
    ylist=[r'$%5.3f^{\circ}$'%i for i in ylist]
#    ax.set_yticklabels(ylist)
    ax.tick_params(labelsize=24)
    ax.tick_params(axis='y',labelright=True,labelleft=False,right=True,left=True)
    ax.tick_params(axis='x',top=True)#,left=True)

    ax.set_xticklabels(xlist,{'fontsize':20})
    ax.set_yticklabels(ylist,{'fontsize':20})
    #--- frame
#    DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
    DrawFrame(ax, (0.1,0.16),(0.09,0.03),0.1)

fig.savefig('%s/spatialMap.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)

AttributeErrorTraceback (most recent call last)
<ipython-input-77-e688aaeedc7d> in <module>()
    111 dd={'ROOT':1505,
    112     'exponent':4.0*0.8,
--> 113      'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
    114       swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
    115     'head_width':.15,'head_length':.3,

AttributeError: class tree has no attribute 'pairs'

spatial map associated with the trig. cascade

In [78]:
#--- plot clusters


def Inside(t,(tmin,tmax)):
    if tmin<= t<tmax:
        return True
    return False    

def GetTrigList( nij_trig ):
    d_trig = {}
    for items in nij_trig.itertuples():
        triggerID = items.Parent_id
        d_trig.setdefault( triggerID, [] ).append( items.Event_id )
    return d_trig

    
def PlotArrowsSpace(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
    xlist=[]
    ylist=[]
    mlist=[]
    tlist=[]
    xm=[]
    ym=[]
    m_mother=[]
    tm=[]
    for triggerID in d_trig: #--- mother events
        
        t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
        ( x0, y0 ) = ( swarm_lohi['longitude'].iloc[triggerID], 
                       swarm_lohi['latitude'].iloc[triggerID] )
        m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
        if not ( #Inside(t0,(tmin,tmax) ) and #--- within the given interval
                  Inside(x0,(xmin,xmax)) and 
                  Inside(y0,(ymin,ymax)) ):
             continue

        for daughter_id in d_trig[triggerID]: #--- children
            t1 = df_complete['date'].iloc[ daughter_id ] #--- time
            ( x1, y1 ) = ( swarm_lohi['longitude'].iloc[ daughter_id ], 
                           swarm_lohi['latitude'].iloc[ daughter_id ] )
            m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
                        
            if not ( #Inside(t1,(tmin,tmax) ) and #--- within the given interval
                      Inside(x1,(xmin,xmax)) and 
                      Inside(y1,(ymin,ymax)) ):
                 continue

            if daughter_id in tree.pairs:
                xlist.append( x1 )
                ylist.append( y1 )
                mlist.append( m1 )
                tlist.append( t1 )
                ax.annotate("", (x1,y1), xytext=(x0, y0), zorder=1,
                         textcoords='data', annotation_clip=False, 
#                         arrowprops=dict(width=.2,headlength=.4,shrink=0.5,color='b'))
                            # backgroundcolor="w",
                        arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
                                        shrinkA=shrinkA,
                                        color="b",linewidth="0.3")) 
 
        #--- plot mother
        if triggerID == tree.pairs[0]: # in roott:
            xm.append(x0)
            ym.append(y0)
            m_mother.append(m0)
            tm.append(t0)
        

    #--- plot circles
    df=pd.DataFrame({'x(km)':xlist,'y(km)':ylist,'m':mlist,'t':tlist})
    ax.scatter(df['x(km)'], df['y(km)'],
                s=exponent**(df['m']),
                alpha=.5, zorder=2,# c = df['t']-tm[0],cmap='jet')#,
                facecolors='red',color='black')
    

    for x0,y0,m0 in zip(xm,ym,m_mother):
        plt.scatter(x0,y0,
                    s=exponent**m0,
                    marker='*',zorder=3,
                    facecolors='yellow',color='black',
                    alpha=1.0)
    
# #--------------------------------------
# #----- key: event value: aftershock id
# #--------------------------------------
# d_trig = GetTrigList( nij_trig )

# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# swarm_lohi.reset_index(inplace=True,drop=True)

# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )


# dd={'ROOT':407,
#     'exponent':4.0*2.7,
#     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#     swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}
# dd={'ROOT':325,
#     'exponent':4.0*2.7,
#      'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}

#--- plot arrows
#for ii in xrange(1):
swarm_lohi = GetCartesian( swarm_lohi )
dff = swarm_lohi.sort_values('magnitude',ascending=False)
roott = [i[1] for i in sl]

    
xxx=np.array(sl)

sr=[[j,i] for i,j in zip(root,S)]
sr.sort(reverse=True)
#sr[:10]

#for ROOT in root[cond]:
for ii in sr[:5]:
    ROOT = ii[1] 
#    if not ROOT == 1:
#        continue
    cls_id = xxx[xxx[:,1]==ROOT][0][2]
    tree.pairs = []
    stree[cls_id].plot(ax)
    tree.pairs = [ROOT]+tree.pairs
    
    dd={'ROOT':ROOT,
    'exponent':4.0*0.8,
#     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
    'dx':0.3*max(swarm_lohi['longitude'].max()-swarm_lohi['longitude'].min(),
    swarm_lohi['latitude'].max()-swarm_lohi['latitude'].min()),
    'head_width':.15,'head_length':.3,
     'shrinkA':0}
    exponent=dd['exponent']
    dx = dd['dx']
    head_width = dd['head_width']
    head_length=dd['head_length']
    shrinkA=dd['shrinkA']

#for ii in xrange(1):
#ROOT = dd['ROOT']
    #for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
    item = dff.loc[ROOT]
    if 1:
#    for item in dff.itertuples(): #--- plot clusters associated with large events
#        if item.Index != ROOT:
#            continue
        
        (xc,yc)=(item['longitude'],item['latitude']) 
#        (xc,yc)=(swarm_lohi['longitude'].mean(),swarm_lohi['latitude'].mean())
        tc = item['date'] #pd.to_datetime("2016-09-03 12:02:44.520")
        dt = 1000   #--- setup
    
        
        fig = plt.figure(figsize=(8,8)) #,frameon=False) #, dpi=300)
        ax = fig.add_subplot(111)
        
        #--- plot all events
        ax.scatter(swarm_lohi['longitude'], swarm_lohi['latitude'], # 'x(km)'],swarm_lohi['y(km)'],
                    s=exponent**swarm_lohi['magnitude'],
                    alpha=0.1)
    
        tt = tc #swarm_lohi['date'].loc[key_max]
        xx={0:xc,1:yc}
        PlotArrowsSpace( swarm_lohi, d_trig,
        ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)), 
               (xx[0]-dx,xx[0]+dx),
               (xx[1]-dx,xx[1]+dx) )
    
        
        ax.set_xlim(xx[0]-dx,xx[0]+dx)
        ax.set_ylim(xx[1]-dx,xx[1]+dx)
        
        
        #--- set tick labels
        xlist=np.linspace(xx[0]-dx,xx[0]+dx,4,endpoint=True)
        ax.set_xticks(xlist)
        xlist=[r'$%5.3f^{\circ}$'%i for i in xlist]
        #
        ylist=np.linspace(xx[1]-dx,xx[1]+dx,4,endpoint=True)
        ax.set_yticks(ylist)
        ylist=[r'$%5.3f^{\circ}$'%i for i in ylist]
    #    ax.set_yticklabels(ylist)
        ax.tick_params(labelsize=24)
        ax.tick_params(axis='y',labelright=True,labelleft=False,right=True,left=True)
        ax.tick_params(axis='x',top=True)#,left=True)
    
        ax.set_xticklabels(xlist,{'fontsize':20})
        ax.set_yticklabels(ylist,{'fontsize':20})
        #--- frame
    #    DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
        DrawFrame(ax, (0.1,0.16),(0.09,0.03),0.1)

        fig.savefig('%s/spatialMap.%s.png'%(DIR_OUTPT_figs,ROOT),bbox_inches='tight',dpi=150)

plot separate clusters

In [79]:
#--- plot clusters


def Inside(t,(tmin,tmax)):
    if tmin<= t<tmax:
        return True
    return False    

def GetTrigList( nij_trig ):
    d_trig = {}
    for items in nij_trig.itertuples():
        triggerID = items.Parent_id
        d_trig.setdefault( triggerID, [] ).append( items.Event_id )
    return d_trig

    
def PlotArrowsSpace(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
    xlist=[]
    ylist=[]
    mlist=[]
    tlist=[]
    xm=[]
    ym=[]
    m_mother=[]
    tm=[]
    for triggerID in d_trig: #--- mother events
        
        t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
#        ( x0, y0 ) = ( swarm_lohi['longitude'].iloc[triggerID], 
#                       swarm_lohi['latitude'].iloc[triggerID] )
        ( X0, Y0 ) = ( swarm_lohi['x(km)'].iloc[triggerID], 
                       swarm_lohi['y(km)'].iloc[triggerID] )
        m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
        
#        if not ( #Inside(t0,(tmin,tmax) ) and #--- within the given interval
#                  Inside(X0,(xmin,xmax)) and 
#                  Inside(Y0,(ymin,ymax)) ):
#             continue

        for daughter_id in d_trig[triggerID]: #--- children
            t1 = df_complete['date'].iloc[ daughter_id ] #--- time
#            ( x1, y1 ) = ( swarm_lohi['longitude'].iloc[ daughter_id ], 
#                           swarm_lohi['latitude'].iloc[ daughter_id ] )
            ( X1, Y1 ) = ( swarm_lohi['x(km)'].iloc[ daughter_id ], 
                           swarm_lohi['y(km)'].iloc[ daughter_id ] )
            m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
                        
#            if not ( #Inside(t1,(tmin,tmax) ) and #--- within the given interval
#                      Inside(X1,(xmin,xmax)) and 
#                      Inside(Y1,(ymin,ymax)) ):
#                 continue

            if daughter_id in d_trig[ROOT]:
                xlist.append( X1 )
                ylist.append( Y1 )
                mlist.append( m1 )
                tlist.append( t1 )
                ax.annotate("", (X1,Y1), xytext=(X0, Y0), zorder=1,
                         textcoords='data', annotation_clip=False, 
#                         arrowprops=dict(width=.2,headlength=.4,shrink=0.5,color='b'))
                            # backgroundcolor="w",
                        arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
                                        shrinkA=shrinkA,
                                        color="b",linewidth="0.3")) 
 
        #--- plot mother
        if triggerID == ROOT: #tree.pairs[0]: # in roott:
            xm.append(X0)
            ym.append(Y0)
            m_mother.append(m0)
            tm.append(t0)
        

    #--- plot circles
    df=pd.DataFrame({'x(km)':xlist,'y(km)':ylist,'m':mlist,'t':tlist})
    ax.scatter(df['x(km)'], df['y(km)'],
                s=exponent**(df['m']),
                alpha=.5, zorder=2,# c = df['t']-tm[0],cmap='jet')#,
                facecolors='red',color='black')

    for x0,y0,m0 in zip(xm,ym,m_mother):
        plt.scatter(x0,y0,
                    s=exponent**m0,
                    marker='*',zorder=3,
                    facecolors='yellow',color='black',
                    alpha=1.0)
    print m_mother[0]
    return (df['x(km)']-xm[0] + 1j*(df['y(km)']-ym[0])).abs(),df['t']-tm[0]
# #--------------------------------------
# #----- key: event value: aftershock id
# #--------------------------------------
d_trig = GetTrigList( nij_trig )

# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
#                              'magnitude', 
#                              ( mc, sys.maxint ) ) 
# swarm_lohi.reset_index(inplace=True,drop=True)

#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )


# dd={'ROOT':407,
#     'exponent':4.0*2.7,
#     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#     swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}
# dd={'ROOT':325,
#     'exponent':4.0*2.7,
#      'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
#     'head_width':.15,'head_length':.3,
#      'shrinkA':0}

#--- plot arrows
#for ii in xrange(1):
swarm_lohi = GetCartesian( swarm_lohi )
#dff = swarm_lohi.sort_values('magnitude',ascending=False)
#roott = [i[1] for i in sl]

    
#xxx=np.array(sl)

#sr=[[j,i] for i,j in zip(root,S)]
#sr.sort(reverse=True)
#sr[:10]

#for ROOT in root[cond]:
if 1: #for ii in sr[:20]:
#     ROOT = ii[1] 
#     if not ROOT == 1:
#         continue
#     cls_id = xxx[xxx[:,1]==ROOT][0][2]
#     tree.pairs = []
#     stree[cls_id].plot(ax)
#     tree.pairs = [ROOT]+tree.pairs
    ROOT = 5 #1509 #5
    dd={'ROOT':ROOT,
    'exponent':4.0*0.8,
#     'dx':max(swarm_lohi.iloc[tree.pairs]['longitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['longitude'],
#      swarm_lohi.iloc[tree.pairs]['latitude'].max()-swarm_lohi.iloc[tree.pairs[0]]['latitude']),
    'dx':0.2*max(swarm_lohi['x(km)'].max()-swarm_lohi['x(km)'].min(),
    swarm_lohi['y(km)'].max()-swarm_lohi['y(km)'].min()),
    'head_width':.15,'head_length':.3,
     'shrinkA':0}
    exponent=dd['exponent']
    dx = dd['dx']
    head_width = dd['head_width']
    head_length=dd['head_length']
    shrinkA=dd['shrinkA']

#for ii in xrange(1):
#ROOT = dd['ROOT']
    #for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
    item = swarm_lohi.loc[ROOT]
    if 1:
#    for item in dff.itertuples(): #--- plot clusters associated with large events
#        if item.Index != ROOT:
#            continue
        
        (xc,yc)=(item['x(km)'],item['y(km)']) 
#        (xc,yc)=(swarm_lohi['longitude'].mean(),swarm_lohi['latitude'].mean())
        tc = item['date'] #pd.to_datetime("2016-09-03 12:02:44.520")
        dt = 1000   #--- setup
    
        
        fig = plt.figure(figsize=(8,8)) #,frameon=False) #, dpi=300)
        ax = fig.add_subplot(111)
        
        #--- plot all events
        ax.scatter(swarm_lohi['x(km)'], swarm_lohi['y(km)'], # 'x(km)'],swarm_lohi['y(km)'],
                    s=exponent**swarm_lohi['magnitude'],
                    alpha=0.1)
    
        tt = tc #swarm_lohi['date'].loc[key_max]
        xx={0:xc,1:yc}
        junk,junk_t=PlotArrowsSpace( swarm_lohi, d_trig,
        ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)), 
               (xx[0]-dx,xx[0]+dx),
               (xx[1]-dx,xx[1]+dx) )
    
        
        ax.set_xlim(xx[0]-dx,xx[0]+dx)
        ax.set_ylim(xx[1]-dx,xx[1]+dx)
        
        
        #--- set tick labels
        xlist=np.linspace(xx[0]-dx,xx[0]+dx,4,endpoint=True)
        ax.set_xticks(xlist)
        xlist=[r'$%s$'%(int(i)) for i in xlist]
        #
        ylist=np.linspace(xx[1]-dx,xx[1]+dx,4,endpoint=True)
        ax.set_yticks(ylist)
        ylist=[r'$%s$'%(int(i)) for i in ylist]
    #    ax.set_yticklabels(ylist)
        ax.tick_params(labelsize=24)
        ax.tick_params(axis='y',labelright=True,labelleft=False,right=True,left=True)
        ax.tick_params(axis='x',top=True)#,left=True)
    
        ax.set_xticklabels(xlist,{'fontsize':20})
        ax.set_yticklabels(ylist,{'fontsize':20})
        #--- frame
    #    DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
        DrawFrame(ax, (0.1,0.16),(0.09,0.03),0.1)

fig.savefig('%s/spatialMap.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)

KeyErrorTraceback (most recent call last)
<ipython-input-79-e40c119402e6> in <module>()
    184         ( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)),
    185                (xx[0]-dx,xx[0]+dx),
--> 186                (xx[1]-dx,xx[1]+dx) )
    187 
    188 

<ipython-input-79-e40c119402e6> in PlotArrowsSpace(df_complete, d_trig, (tmin, tmax), (xmin, xmax), (ymin, ymax))
     51 #                 continue
     52 
---> 53             if daughter_id in d_trig[ROOT]:
     54                 xlist.append( X1 )
     55                 ylist.append( Y1 )

KeyError: 5
In [80]:
ROOT = 1505 #1

fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlim(100,1e6)
ax.set_ylim(1e-1,1e3)
ax.set_xscale('log')
ax.set_yscale('log')
    
for keys in stree:
    tree_obj = stree[ keys ]
    if tree_obj.data != ROOT:
        continue
    tree.rtm_dic = {} #--- initialize distance-time list
    tree_obj.get_coords() #--- update
    
    t_root = tree_obj.coord[ 0 ] #--- relative to the root
    r_root = np.array(tree_obj.coord[ 1 ][:2] )
    delta_t = []
    delta_r = []
    mmm=[]
    for key in tree.rtm_dic:
        t_afs = tree.rtm_dic[key][0]
        r_afs = np.array(tree.rtm_dic[key][1][:2])
#         if key == 5:
#             print tree.rtm_dic[key][0]
#             print tree.rtm_dic[key][1][:2]
#             print tree.rtm_dic[key][2]
        delta_t += [ ( t_afs - t_root ).total_seconds() ]
        delta_r += [ ( ( ( r_afs - r_root ) ** 2 ).sum() ) ** 0.5 ]
        mmm += [tree.rtm_dic[key][ 2 ]]

ax.scatter(delta_t,delta_r,s=3**np.array(mmm),alpha=0.1)
        

NameErrorTraceback (most recent call last)
<ipython-input-80-6eb46335e5fd> in <module>()
     31         mmm += [tree.rtm_dic[key][ 2 ]]
     32 
---> 33 ax.scatter(delta_t,delta_r,s=3**np.array(mmm),alpha=0.1)
     34 

NameError: name 'delta_t' is not defined
In [81]:
#print tree.rtm_dic

# swarm_lohi.sort_values(by='date',inplace=True)            
# coord_length=swarm_lohi['longitude'].apply(lambda x:[x])+swarm_lohi['latitude'].apply(lambda x:[x])
# Magnitude = np.exp(swarm_lohi['magnitude'])/20000

KEYS=tree.rtm_dic.keys()
KEYS.sort()

coordd={}
Magnitudee={}
for i in KEYS:
    coordd[i]=tree.rtm_dic[i][1][:2]
    Magnitudee[i]=np.exp(tree.rtm_dic[i][2])/1000
makeDATA( open( '%s/map2.xyz'%DIR_OUTPT, 'w' ), None ).Map( coordd, 'ITIME=%s'%0, Magnitudee) #--- instantanuous
EVENT_ID = Magnitudee.keys()
EVENT_ID = map(int,EVENT_ID)
EVENT_ID.sort()
cords={}
mag={}
t0 = None
os.system('rm %s/temporalMap2.xyz'%DIR_OUTPT)
for j_event in EVENT_ID:
    t = swarm_lohi['date'].iloc[ j_event ]
    if not t0:
        t0 = t
    assert t >= t0, 't0=%s,t=%s'%( t0, t )
    #--- temporal map 
    cords[ j_event ] = coordd[ j_event ]
    mag[ j_event ] = Magnitudee[ j_event ]
    itime=swarm_lohi[ 'date' ].iloc[ j_event ]
    tstr='%s-%s'%(itime.strftime("%x"), itime.strftime("%X"))
    makeDATA( open( '%s/temporalMap2.xyz'%DIR_OUTPT, 'a' ), None ).Map( cords, 'ITIME=%s'%itime, mag ) #--- instantanuous
    t0 = t

NameErrorTraceback (most recent call last)
<ipython-input-81-0066fe586382> in <module>()
     13     coordd[i]=tree.rtm_dic[i][1][:2]
     14     Magnitudee[i]=np.exp(tree.rtm_dic[i][2])/1000
---> 15 makeDATA( open( '%s/map2.xyz'%DIR_OUTPT, 'w' ), None ).Map( coordd, 'ITIME=%s'%0, Magnitudee) #--- instantanuous
     16 EVENT_ID = Magnitudee.keys()
     17 EVENT_ID = map(int,EVENT_ID)

NameError: name 'makeDATA' is not defined
In [82]:
swarm_loh = swarm_lohi.iloc[6:1510]
r0=np.array(swarm_lohi.iloc[5][['x(km)','y(km)']])
t0 = swarm_lohi.iloc[5][['date']]
m0 = swarm_lohi.iloc[5][['magnitude']]
dt=swarm_loh['date'].apply(lambda x:(x-t0))
dr=swarm_loh['x(km)'].apply(lambda x:x-r0[0])+swarm_loh['y(km)'].apply(lambda x:1j*(x-r0[1]))
dr=dr.abs()
dt_sec = dt['date'].apply(lambda x:x.total_seconds()/3600)
print m0
magnitude    3
Name: 5, dtype: object
In [83]:
fig=plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xlim(1e-1,1e2)
ax.set_ylim(1e-1,1e2)
ax.set_xscale('log')
ax.set_yscale('log')
ax.scatter(dt_sec,dr,s=3.0**(swarm_loh['magnitude']),alpha=0.2)
Out[83]:
<matplotlib.collections.PathCollection at 0x7fe7365886d0>
In [84]:
#--- r vs t (including indirect aftershocks)

dd = {}    
for cls_id in stree:
    index_m = np.array(m).searchsorted(stree[cls_id].coord[2]) - 1 #--- group clusters based on magnitude range
    dd.setdefault(index_m,[]).append(cls_id)

colors = ['black','red','green','blue','cyan','brown']

for group_m, col in zip(dd, colors ): #--- loop over magnitude ranges
    #--- plot
    fig = plt.figure(figsize=(4,4))
    ax = fig.add_subplot(111)
    ax.set_xlim(1,1e5)
    ax.set_ylim(1e-3,1e2)
    ax.set_xscale('log')
    ax.set_yscale('log')


    cls_ids = dd[ group_m ] #--- cluster ids for a given range
    delta_t = []
    delta_r = []
    for cls_id in cls_ids: #--- loop over each cluster
        tree.rtm_dic = {} #--- initialize distance-time list
        stree[cls_id].get_coords() #--- update
        t_root = stree[cls_id].coord[ 0 ] #--- relative to the root
        r_root = np.array( stree[cls_id].coord[ 1 ] )
        for key in tree.rtm_dic:
            t_afs = tree.rtm_dic[key][0]
            r_afs = np.array(tree.rtm_dic[key][1])
            delta_t += [ ( t_afs - t_root ).total_seconds() ]
            delta_r += [ ( ( ( r_afs - r_root ) ** 2 ).sum() ) ** 0.5 ]
        
    
    ax.scatter(delta_t,delta_r,color=col,alpha=0.05)

plt.show()
In [85]:
t0
Out[85]:
date    2019-07-04 17:35:52
Name: 5, dtype: object
In [86]:
#--- r vs t (including indirect aftershocks)

dd = {}    
for cls_id in stree:
    index_m = np.array(m).searchsorted(stree[cls_id].coord[2]) - 1 #--- group clusters based on magnitude range
    dd.setdefault(index_m,[]).append(cls_id)

colors = ['black','red','green','blue','cyan','brown']

for group_m, col in zip(dd, colors ): #--- loop over magnitude ranges
    #--- plot
    fig = plt.figure(figsize=(4,4))
    ax = fig.add_subplot(111)
    ax.set_xlim(1,1e5)
    ax.set_ylim(1e-3,1e2)
    ax.set_xscale('log')
    ax.set_yscale('log')


    cls_ids = dd[ group_m ] #--- cluster ids for a given range
    delta_t = []
    delta_r = []
    for cls_id in cls_ids: #--- loop over each cluster
        tree.rtm_dic = {} #--- initialize distance-time list
        stree[cls_id].get_coords() #--- update
        t_root = stree[cls_id].coord[ 0 ] #--- relative to the root
        r_root = np.array( stree[cls_id].coord[ 1 ] )
        for key in tree.rtm_dic:
            t_afs = tree.rtm_dic[key][0]
            r_afs = np.array(tree.rtm_dic[key][1])
            delta_t += [ ( t_afs - t_root ).total_seconds() ]
            delta_r += [ ( ( ( r_afs - r_root ) ** 2 ).sum() ) ** 0.5 ]
        
    
    ax.scatter(delta_t,delta_r,color=col,alpha=0.05)

plt.show()
In [ ]: